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Abstract 

A  three  dimensional,  second  order  finite  difference  method  was  used  to  create  synthetic 
seismograms  for  elastic  wave  propagation  in  heterogeneous  media.  These  synthetic  seismograms  are 
used  to  model  rough  seafloor,  the  shallow  crust,  or  complex  structural  and  stratigraphic  settings  with 
strong  lateral  heterogeneities.  The  finite  difference  method  is  preferred  because  it  allows  models  of 
any  complexity  to  be  generated  and  includes  all  multiple  scattering,  wave  conversion  and  diffraction 
effects.  The  method  uses  a  fully  staggered  grid  as  developed  by  Virieux  (1986).  Wavefront  snapshots 
and  time  series  output  allow  the  scattering  and  focussing  of  different  wave  modes  with  direction  to 
be  visualized. 

The  extensive  calculations  required  for  realistic  size  models  stretches  the  resources  of  serial 
computers  like  the  VAX  8800.  On  the  Connection  Machine,  a  massively  parallel  computer,  the 
finite  difference  grid  can  be  directly  mapped  onto  the  virtual  processors,  reducing  the  nested  time 
and  space  loops  in  the  serial  code  to  a  single  time  loop.  As  a  result,  the  computation  time  is  reduced 
dramatically.  .  _ 

■  -  >  -  '  '■  •*-'  /-  'Y'-  * 


IV 


1  Introduction 


Most  computers  used  today  are  based  on  the  von  Neuman  (serial)  architecture.  Serial  computers 
have  a  single  central  processing  unit  (CPU)  and  memory;  data  from  memory  is  passed  to  the  CPU 
one  data  element  at  a  time  and  instructions  are  executed  in  sequence.  These  computers  are  termed 
single-instruction,  single-datastream  (SISD)  architecture  since  one  instruction  at  a  time  is  executed 
on  a  single  data  item.  Although  the  efficiency  of  SISD  computers  may  be  enhanced  by  hardware 
techniques  such  as  pipelining  and  vector  processors,  the  time  required  to  complete  computationally 
intensive  programs  is  often  prohibitive. 

Parallel  systems  can  be  grouped  into  two  broad  categories:  single-instruction,  multiple-datastream 
(SIMD)  and  multiple-instruction,  multiple-datastream  (MIMD).  SIMD  computers  are  characterized 
by  a  large  number  of  simple  processors,  each  with  its  own  local  memory  connected  by  a  network  com¬ 
munications  system.  MIMD  computers  typically  have  a  smaller  number  of  conventional  processors 
with  shared  memory  which  execute  portions  of  a  program  concurrently. 

In  an  application  such  as  synthetic  seismograms,  which  requires  both  large  volumes  of  data 
and  extensive  calculations,  the  time  required  to  cycle  all  the  data  through  one  CPU  places  severe 
limitations  on  the  size  of  models  which  can  be  run  on  a  serial  computer.  This  report  documents  the 
procedure  used  to  implement  a  software  system  to  compute  3-dimensional  synthetic  seismograms 
by  the  finite  differences  method  on  the  Connection  Machine,  a  SIMD  computer.  The  advantage  of 
parallel  processing  on  a  SIMD  computer  is  that  each  node  of  the  3D  model  can  be  assigned  to  a 
single  processor.  Calculations  are  then  performed  simultaneously  on  all  grid  points.  An  operation 
which  would  normally  be  performed  within  a  repetitive  loop  is  replaced  by  a  single  operation  on 
many  processors  acting  in  parallel. 

Two  dimensional  synthetic  seismogram  models  were  previously  coded  in  FORTRAN  77  to  run 
on  VAX  11/780,  VAX  8800,  Cyber  205  and  Cray  XMP-12  computers  (Hunt  and  Stephen,  1986).  A 
FORTRAN  77  program  for  a  three-dimensional  solution  to  the  elastic  wave  equation  by  the  method 
of  finite  differences,  which  runs  on  the  VAX  8800,  was  used  as  the  basis  for  the  parallel  programs 
written  in  C*.  The  parallel  program  is  called  wave_3d. 


1.1  The  Connection  Machine 

The  Connection  Machine  (CM)  is  a  massively  parallel  computer  developed  by  Thinking  Ma¬ 
chines,  in  Cambridge,  MA.  Massively  parallel,  or  fine-grain  computers  have  many  simple  processors, 
each  with  its  own  local  memory.  The  Connection  Machines  used  in  this  project  are  located  at 
the  Naval  Research  Laboratory  (NRL)  Connection  Machine  Facility  and  at  the  Northeast  Paral¬ 
lel  Architectures  Center  (NPAC)  at  Syracuse  University.  NRL  has  two  CM-2  machines,  one  with 
8k  processors  (Bambi)  and  one  with  16k  processors  (Godzilla).  NPAC  also  has  two  Connection 
Machines,  a  32k  CM-1  (CUBE)  and  a  32k  CM-2  (SON-OF-CUBE). 

The  Connection  Machine  processors  are  divided  into  groups  of  8k  or  16k  processors.  A  single 
Connection  Machine  can  have  one,  two  or  four  groups  of  processors.  The  16k  NRL  CM-2  has  two 
groups  of  8k  processors;  the  NPAC  32k  CM-2  has  4  groups  of  8k  processors.  Each  group  of  processors 
is  associated  with  a  sequencer  which  interprets  the  instructions  sent  by  a  front-end  computer.  In 
order  to  use  the  Connection  Machine  system,  the  front-end  must  logically  attach  itself  to  one  or 
more  sequencers. 
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Each  processor  on  the  Connection  Machine  model  CM-2  has  64K  bits  (2048  words)  of  memory; 
each  CM-1  processor  has  4k  bits  of  memory.  The  Connection  Machine  processors  are  connected  by 
a  network  communications  system  so  that  any  processor  can  communicate  with  any  other  processor 
via  a  routing  device  wired  in  an  n-cube  pattern  (Hillis,1987).  The  CM  can  be  configured  in  software 
for  virtual  to  physical  processor  ratios  that  are  a  power  of  2,  as  long  as  there  is  enough  memory  per 
virtual  processor  to  accomodate  the  required  calculations  and  variables. 

The  memory  requirements  of  the  wave_3d  program  limits  the  maximum  virtual  processor  ratio 
to  16:1.  On  the  16k  NRL  Connection  Machine  (Godzilla),  the  maximum  number  of  virtual  processors 
available  for  the  3D  grid  is  262,144  (64x64x64  or  256k);  on  the  32k  NPAC  CM-2  Connection  Machine 
(SON-OF-CUBE),  the  maximum  number  of  grid  points  is  524,288  (64x64x128  or  512k). 


1.2  Programming  the  Connection  Machine 


The  CM  is  connected  by  a  high  speed  bus  to  a  conventional  serial  computer  which  serves  as  the 
user-interface.  Program  development  is  accomplished  on  the  front-end  computer  using  the  editors 
and  compilers  of  the  front-end  processor.  All  of  the  Connection  Machines  used  in  this  project  support 
C*,  C/Paris,  *Lisp  and  CM  Fortran. 


The  NRL  Connection  Machines  Bambi  (8k)  and  Godzilla  (16k)  are  connected  to  four  front-end 
computers: 


1.  cmvax 

2.  cmsun 

3.  THINK75 

4.  THINK40 


(VAX  8800) 
(SUN  4/280) 
(Symbolics  3675) 
(Symbolics  3640) 


The  NRL  CMs  are  interfaced  to  the  front-ends  as  follows: 


Interface  cmvax  cmsun  THINK75  THINK40 

0  Bambi  Bambi  Godzilla  Bambi 

1  Godzilla  Godzilla 


The  VAX  front-end  (cmvax)  runs  the  ULTRIX  operating  system;  the  SUN  front-end  (cmsun) 
runs  the  UNIX  operating  system. 

The  two  Connection  Machines  at  NPAC  (CUBE  and  SON-OF-CUBE)  are  connected  to  a  VAX 
front-end  (cmx)  which  runs  ULTRIX  2.2,  a  4.2BSD-based  system. 

The  parallel  program  is  downloaded  onto  the  CM  at  runtime.  The  synthetic  seismogram  program, 
wave_3d,  is  written  in  C*,  an  extension  of  the  C  programming  language,  developed  to  allow  parallel 
execution.  The  wave_3d  program  contains  both  serial  and  parallel  code. 


2  The  Model 


The  finite  difference  method  is  used  to  create  synthetic  seismograms  for  models  containing  lateral 
heterogeneity;  this  method  is  preferred  because  it  allows  models  of  any  complexity  to  be  generated 
and  includes  all  multiple  scattering,  wave  conversion,  and  diffraction  effects.  Most  analytical  mod¬ 
elling  techniques  use  assumptions,  such  as  primary  scattering  (Born  and  Rytov  approximations), 
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which  limit  their  applicability  to  weak  heterogeneity  situations.  Modelling  of  a  rough  seafloor,  the 
shallow  crust,  or  complex  structural  and  stratigraphic  settings,  however,  must  be  able  to  handle 
strong  lateral  heterogeneities.  In  addition,  the  finite  difference  technique  can  be  used  for  all  ratios 
of  scatterer  size  to  wavelength. 


2.1  Background 


The  finite  difference  method  involves  the  spatial  and  temporal  discretization  of  the  wave  equation 
on  a  regular  grid  which  represents  the  model  of  interest.  The  method  has  been  successfully  applied 
to  an  extremely  wide  range  of  seismic  wave  propagation  problems  including  earthquake  seismology 
(Alterman  and  Karal,  1968;  Frankel  and  Clayton,  1986;  Toksoz  et  al.,  1988),  marine  refraction 
(Stephen  1983;  Dougherty  and  Stephen,  1988),  reflection  seismology  (Kelly  et  al.  1976;  Virieux, 
1986,  Fornberg,  1987),  VSP  (Stephen,  1984),  and  full  waveform  acoustic  logging  (Stephen  et  al., 
1985).  A  number  of  formulations  have  been  utilized  for  these  applications,  including  second  and 
fourth  order  formulations,  and  the  pseudospectral  method.  Fornberg  (1987)  presents  a  comparison 
of  these  different  methods.  The  higher  order  schemes  are  more  accurate,  require  fewer  grid  points  per 
wavelength,  and  are  more  computationally  demanding.  The  lower  order  schemes  are  more  efficient 
computationally,  but  require  more  grid  points  per  wavelength,  and  therefore  more  computer  memory. 
The  second  order  scheme  has  been  fully  validated  against  analytic  methods  (Stephen,  1983,  1988) 
and  can  handle  complicated  interfaces  and  boundaries.  In  addition,  the  use  of  fully  staggered  grids 
for  displacements  and  stresses  (Virieux,  1986;  Dougherty  and  Stephen,  1988)  improves  the  accuracy 
and  stability  of  the  formulation  at  no  additional  cost  in  computation  or  memory.  The  second  order 
scheme  has  been  successfully  applied  to  2-D  heterogeneous  media  by  Dougherty  and  Stephen  (1988), 
and  3-D  media  by  Etgen  and  Yomagida  (1988).  Because  the  primary  interest  is  in  modelling  rough 
interface  effects,  even  the  pseudospectral  method  will  need  fine  grid  spacing  to  accurately  represent 
such  interfaces.  Therefore,  the  memory  requirement  differences  between  the  lower  order  methods 
and  the  pseudospectral  method  we  not  as  great  as  for  other  types  of  applications. 


2.2  Approach 

The  system  of  equations  we  wish  to  solve  is  given  by: 
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Where  u,v,w  represent  the  x,y,z  displacements,  X  and  p  are  the  Lame  parameters,  p  is  the 
density,  r jj  the  stress  terms,  and  t  represents  time.  For  heterogeneous  media,  the  Lame  parameters 
and  density  all  vary  spatially  in  three  dimensions.  A  unique  staggered  spatial  grid  can  be  constructed 
to  solve  this  system  with  centered  differences.  This  grid  is  a  direct  extension  of  the  two  dimensional 
grid  given  in  Virieux  (1986). 

In  order  to  avoid  grid  dispersion  in  the  model,  the  highest  frequency  component  to  be  modelled 
must  be  sampled  by  10  -  30  grid  points  per  wavelength  (Kelly  et  al.,  1976;  Stephen,  1983). 

A  dilatational  or  explosive  point  source  is  introduced  into  the  grid  using  the  method  developed 
by  Nicoletis  (1981).  This  method  computes  a  force  distribution  over  a  finite  volume  of  the  three 
dimensional  grid  to  represent  a  dilatational  source.  The  accuracy  of  this  representation  increases 
as  the  grid  volume  on  which  it  is  imposed  is  increased.  The  time  dependency  of  the  source  takes 
the  form  of  the  derivative  of  the  Gaussian  distribution  with  any  given  center  frequency  (Kelly  et  al. 
1976;  Stephen  et  al.,  1985). 

Boundary  conditions  for  this  3-dimensional  staggered  grid  formulation  are  fairly  straight  forward. 
The  source  is  introduced  inside  the  model,  and  the  six  planes  which  define  the  model  limits  are 
treated  as  absorbing  boundaries  using  the  telegraph  equation  in  a  zone  around  each  of  the  planes 
(Levander,  1985;  Dougherty  and  Stephen,  1988).  Because  the  finite  difference  model  is  formulated 
for  heterogeneous  media,  interfaces  within  the  model  do  not  require  any  specific  boundary  conditions 
to  be  coded.  Interfaces  of  any  complexity  are  handled  implicitly  (with  the  caveat  that  steeply  dipping 
interfaces  must  be  adequately  sampled  by  the  spatial  grid). 

The  memory  and  computational  requirements  for  numerical  modelling  in  three  dimensions  are 
great.  For  the  second  order  three-dimensional  code  which  we  have  developed,  twelve  (12)  real 
variables  must  be  stored  for  each  grid  point  in  the  heterogeneous  section  of  the  model  (three  dis¬ 
placement  values  at  three  time  steps,  plus  two  elastic  constants  and  density).  In  order  to  reduce  the 
memory  requirements  on  the  VAX,  the  model  assumes  a  heterogeneous  zone  sandwiched  between 
two  homogeneous  zones.  As  an  example  of  the  memory  and  computational  requirements  on  a  serial 
machine,  a  model  containing  80  x  80  x  50  grid  points  needs  12.5  Mbytes  of  central  memory  and  the 
run  time  for  400  time  steps  was  about  20  hours  on  a  VAX  8800.  These  requirements  can  be  reduced 
by  reducing  the  heterogeneous  zone  to  as  small  a  volume  as  needed  for  the  models  of  interest  or  by 
segmenting  the  model  into  several  runs. 

Fortunately,  computer  hardware  advances  in  parallel  processing  can  solve  this  problem.  On  a 
massively  parallel  processing  machine  the  finite  difference  grid  can  be  directly  mapped  onto  the 
processors,  reducing  the  nested  time  and  space  loops  in  the  serial  code  to  a  single  time  loop  in  the 
parallel  code.  As  a  result,  the  computation  time  is  reduced  dramatically.  Other  parallel  machines  are 
designed  with  fewer  more  powerful  processors.  On  these  types  of  machines,  each  processor  might  run 
the  finite  difference  code  for  a  segment  of  the  model,  and  communicate  displacement  values  at  the 
end  of  each  time  step.  Any  parallel  computer  can  greatly  reduce  computation  time  and,  therefore, 
make  three  dimensional  numerical  modelling  feasible. 
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3  Running  the  synthetic  seismogram  program 


This  section  describes  the  steps  necessary  to  define  the  model,  run  the  synthetic  seismogram 
program  and  view  the  three-dimensional  model  created  by  the  program  wave -3d.  Although  some 
familiarity  with  the  UNIX  operating  system  is  helpful,  the  following  sections  provide  a  complete 
description  of  the  steps  required. 

To  run  program  wave_3d: 

1.  Log  on  to  the  front-end  via  telnet  (from  the  WHOI  RED  VAX) 

2.  Define  the  model  parameters  in  the  file  xmodel.dat 

3.  Run  the  wave_3d  program  on  the  Connection  Machine 

4.  View  the  results  by  either: 

(a)  copying  the  ASCII  output  files  to  the  VAX  (via  ftp)  and  using  the  SNAPSHOT  program 
on  the  WHOI  THISLE  VAX  Workstation,  or 

(b)  running  wave-3d  with  the  frame  buffer  option  set  and  viewing  the  results  on  the  CM 
Graphic  Display  System  (GDS)  high  resolution  color  monitor  (note  that  there  should  be 
someone  at  NRL  or  NPAC  who  can  physically  watch  the  frame  buffer  monitor!) 

3.1  Logging  on  via  telnet 

Access  to  the  Connection  Machine  front-end  computers  at  NRL  and  NPAC  is  via  the  ARPA 
Internet  link  using  telnet  (figure  1). 

1 .  log  on  to  the  WHOI  Red  VAX 

To  access  the  Connection  Machine  at  NRL: 

RED-$  telnet  134.207.7.12  (for  the  VAX  interface  -  cmvax) 

RED.S  telnet  134.207.7.4  (for  the  SUN  interface  -  cmsun) 

To  access  the  Connection  Machine  at  NPAC: 

RED-$  telnet  cmx.npac.syr.edu  (for  the  VAX  interface  -  cmx) 

Note:  when  using  telnet,  you  could  get  the  following  message: 

REDJ$  telnet  134.207.7.12 

Trying. . .Attempt  to  connect  to  foreign  host  failed:  Hetwork  is  unreachable 
RED-S  <<<  at  this  point,  just  wait  and  try  again  later!  >>> 

2.  log  on  to  cmvax,  cmsun  or  cmx;  after  logging  on,  one  of  the  following  prompts  will  appear: 

cmvax'/.  (if  logged  on  to  the  VAX  front-end  at  NRL) 
cmsun%  (if  logged  on  to  the  SUN  front-end  at  NRL) 
cmx'/.  (if  logged  on  to  the  VAX  front-end  at  NPAC) 
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Figure  la.  Logging  on  to  the  NRL  SUN  front-end  (emsun) 

REDJ*  telnet  134.207.7 .4 

Trying . . .  Open 
Connected  to  134.207.7.4. 

Escape  character  is 

SunOS  UNIX  (emsun) 
login:  alien 

Password:  <<<type  password  here  >>> 

Last  login:  Thu  Aug  31  08:37:32  irom  RED.WH0I.EDU 

SunOS  Release  4.0  (CMSUH)  #64:  Tue  Aug  22  12:43:06  EDT  1989 

URL  Connection  Machine  Facility  -  Sun  4/280 

<<<  system  messages  here  »> 

TERM  =  (vtlOO)  <<<  hit  return  here  >>> 
alien .  cm  sun  7.  to 

Connection  closed  to  134.207.7.4 

R£D_*  «<  back  to  DCL  prompt  on  RED  VAX  >>> 

Figure  lb.  Logging  on  to  the  NRL  VAX  front-end  (cmvax). 

RED-*  telnet  134.207.7.12 

Trying...  Open 
Connected  to  134.207.7.12. 

Escape  character  is  '*]*. 

cmvax  login:  alien 

Password:  <<<  type  password  here  >» 

Last  login:  Tue  Aug  8  11:31:48  irom  128.128.16.79 

Ultrix-32  V3.0  (Rev  64)  System  #7:  Sat  Jul  8  17:34:53  EDT  1989 

HRL  Connection  Machine  Facility  -  VAX  8800 

«<  system  messages  here  »> 

TERM  =  (vtlOO)  «<  hit  return  here  >» 
alien .  cmvax1/.  lo 

Connection  closed  to  134.207.7.12 

RED_t  «<  back  to  DCL  prompt  on  RED  VAX  »> 

Figure  lc.  Logging  on  to  the  NPAC  VAX  front-end  (cmx). 

RED.*  telnet  cmx.npac.syT.edu 

Trying . . . 

Open 

Connected  to  cmx.npac.syr.edu. 

Escape  character  is 


Ultrix-32  V3.0  (Rev  64)  (cmx.npac.syr.edu) 
login:  j alien 

Password:  <<<  type  password  here  >>> 

Last  login:  Tue  Sep  E  09:17:39  from  RED.VH0I.EDU 

Ultrix-32  V3.0  (Rev  64)  System  #1:  Mon  Jul  17  11:24:69  EDT  1989 

«<  system  messages  here  >» 

Tue  Sep  E  09:17:40  EDT  1989 

cmx1/  lo 

Connection  closed  to  cmx.npac.syr.edu 
RED.*  <<<  back  to  DCL  prompt  on  RED  VAX  >>> 


Figure  1:  Example  of  logging  on  and  off  the  NRL  and  NPAC  Connection  Machines 


6 


3.2  Defining  the  model  input  parameters 


Input  to  program  wave_3d  is  via  the  ASCII  file  xmodel.dat.  The  contents  of  xmodel.dat 
are  as  follows: 
nt 

dx,dy,dz,dt 

snap-file, snap-fb 

out-fb,out_xy,out.yz,out_xz 

isx,isy,isz,fpeak 

nzones 

(for  each  zone:) 

nxl  ,ny  1  ,nz  1  ,nxn,nyn,nzn 

cclamb  ,ccmu  ,trho 

where: 

nt 

dx,dy,dz 
dt 

snap-file 
snap-fb 
out-fb 


out-xy 
out.yz 
outjcz 
i8x,isy,isz 
fpeak 
nzones 
nxl,nyl,nzl 
nxn,nyn,nzn 
cclamb 
ccmu 
trho 

Figure  2a  shows  a  sample  input  data  file  for  the  model  illustrated  in  figure  2b.  Note  that  when 
defining  the  boundaries  of  the  model  areas,  the  indices  (nxl,nyl,nzl,nxn,nyn,nzn)  should  range  from 
0  to  63,  rather  than  from  I  to  64. 

Use  the  vi  editor  on  the  front-end  to  edit  this  file  or  use  the  edt  editor  on  the  WHOI  RED  VAX 
and  ftp  the  xmodel.dat  file  to  the  CM  front-end  (see  section  3.4  for  directions  on  transferring  files 
via  ftp). 


=  number  of  time  steps 
=  grid  spacing  in  x,y,z  directions,  in  meters 
=  time  step  increment,  in  seconds 

=  interval  for  output  of  ASCII  files  (number  of  time  steps) 

=  interval  for  output  to  frame  buffer  (number  of  time  steps) 

=  0  for  no  output  to  frame  buffer 

=  1  to  output  xy  plane  to  frame  buffer,  at  snap-fb  intervals 
=  2  to  output  yz  plane  to  frame  buffer,  at  snap-fb  intervals 
=  3  to  output  xz  plane  to  frame  buffer,  at  snap-fb  intervals 
=  1  to  output  xy  divergence  and  curl  z  to  ASCII  files,  at  snap-file  intervals 
=  1  to  output  yz  divergence  and  curl  x  to  ASCII  files,  at  snap-file  intervals 
=  1  to  output  xz  divergence  and  curl  y  to  ASCII  files,  at  snap-file  intervals 
=  location  of  point  source  (grid  points) 

=  peak  frequency  of  source  (hz) 

=  number  of  homogeneous  zones  in  model 
=  start  point  for  zone  in  x,y,z  directions  (meters) 

=  end  point  for  zone  in  x,y,z  directions  (meters) 

=  Lame’s  parameter  -  lambda 
=  Lame’s  parameter  -  mu  (shear  modulus) 

=  density  (kg/m  ) 
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Figure  2a.  Sample  file  xmodel.dat  for  heterogeneous  model. 

This  file  describes  the  parameters  for  the  model  in  figure  2b.  The  3D  grid  is  64x64x64; 
the  point  source  is  located  at  (10,31,10);  there  will  be  no  output  to  the  frame  buffer;  ASCII 
files  created  will  be:  div_xy.dat,  curl_z.dat,  div_xz.dat,  curl_y.dat,  div_yz.dat,  and 
curl-x.dat.  The  model  will  be  calculated  for  250  time  steps  and  data  will  be  output  to  the 
files  at  every  50  time  steps. 


250 

50  50  50  0.007 
50  2 
0  111 
10  31  10  3. 

3 

0  0  0  63  63  31 
2.25e09  0.0  1000. 

0  0  32  31  63  63 
2 . 25e09  0.0  1000. 

32  0  32  63  63  63 
16.0e09  8 . 0e09  2000. 


(number  of  time  steps) 

(dx,  dy,  dz,  dt) 

(snap_files,  snap_fb) 

(out_fb,  out_xy,  out_yz,  out_xz) 

(isx,  isy,  isz,  fpeak) 

(number  ol  zones) 

(nxl,  nyl,  nzl,  nxn,  nyn,  nzn  for  zone  1) 
(cclamb,  ccmu,  trho  for  zone  1) 

(nxl,  nyl,  nzl,  nxn,  nyn,  nzn  f^r  zone  2) 
(cclamb,  ccmu,  trho  for  zone  2) 

(nxl,  nyl,  nzl,  nxn,  nyn,  nzn  for  zone  3) 
(cclamb,  ccmu,  trho  for  zone  3) 


Figure  2b.  Three  dimensional  model  calculated  in  the  examples  presented  in 
this  technical  report. 
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Figure  2:  Sample  input  data  file  (xmodel.dat) 
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3.3  Running  the  program 


Program  wave_3d  is  configured  for  a  model  with  dimensions  64x64x64;  this  model  requires 
262144  (256k)  virtual  processors.  The  program  has  been  optimized  to  run  at  a  maximum  virtual 
processor  ratio  of  16:1.  This  means  that  a  minimum  of  16k  processors  is  necessary  to  run  this  model. 
When  running  the  program  at  NRL,  all  16k  processors  of  the  Godzilla  Connection  Machine  must  be 
allocated.  On  the  32k  NPAC  CM,  only  two  of  the  four  sequencers  on  SON-OF-CUBE  are  needed  to 
run  the  program.  The  program  can  also  be  configured  to  run  a  model  with  dimensions  64x128x64 
on  the  NPAC  32k  CM  if  all  four  sequencers  on  SON-OF-CUBE  are  used  (see  section  4.8). 

To  run  the  wave_3d  program  (figure  3): 

1)  Determine  if  the  required  sequencers  on  the  CM  are  available  using  the  cmfinger  and  cmuser 
commands  (figure  4);  you  will  not  be  allowed  to  attach  to  the  CM  if  it  is  busy 

2a)  on  the  NRL  CM,  type  the  following  command: 

cmvaxV.  cmattach  -b  0.95  -p  16k  -v  256k  -il  wave_3d  >  wave_3d.out  & 

2b)  on  the  NPAC  CM,  use  this  command: 

cmx'/.  cmattach  -b  0.95  -p  16k  -v  256k  -C  S  wave_3d  >  wave_3d.out  & 

The  file  wave.3d.out  will  contain  the  runtime  diagnostic  messages  and  the  timing  information 
after  program  completion.  If  the  CM  is  busy  and  you  wish  to  attach  and  ’wait  for  resources’,  add 
the  -w  flag  to  the  command  line: 

cmvaxtt  cmattach  -b  0.95  -p  16k  -v  256k  -w  -il  wave -3d  >  wave  _3d.  out  & 
or 

cmx'/.  cmattach  -b  0.95  -p  16k  -v  256k  -w  -C  S  wave-3d  >  wave-3d.out  & 

This  option  will  cause  the  cmattach  program  to  wait  (possibly  forever)  for  the  desired  resources 
(namely  both  sequencers  of  Godzilla  or  two  sequencers  of  SON-OF-CUBE)  to  be  freed.  Note  that 
this  is  a  fairly  primitive  implementation  and  simply  causes  the  job  to  be  resubmitted  once  every 
minute.  If  the  sequencers  are  freed  up  during  the  minute  that  the  job  is  waiting,  they  can  be  attached 
by  anyone  else.  If  you  are  in  a  hurry,  it  is  safer,  and  usually  faster,  to  use  cmfinger  to  keep  an  eye 
on  which  sequencers  are  free  and  then  resubmit  the  cmattach  command.  The  c-shell  has  a  history 
command  that  is  useful  for  this: 

cmx%  !cma  (the  most  recent  command  beginning  with  cma  will  be  executed) 

The  output  file  should  have  a  unique  name  on  the  NRL  CM.  If  a  file  with  that  name  exists 
already,  the  following  message  will  appear: 

cnsun'/.  cmattach  -b  0.95  -p  16k  -v  256k  *il  wave-3d  >  wave_3d.out  & 

[1]  13999 

vave.3d.out  tile  exists 

cnsun'/. 

Either  choose  a  different  name  for  the  output  file  on  the  command  line,  or  remove  the  current 
file  and  re-submit  the  attach  command: 
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crnsun'/.  rm  waveJJd.out 
rm:  remove  wave.3d.out?  y 
cmsun'/,  !cma 

On  the  NPAC  CM,  the  file  will  be  overwritten  if  it  already  exists. 

To  have  diagnostic  messages  appear  on  the  terminal  screen  rather  than  in  a  file,  simply  omit  the 
redirection  (>  wave_3d.out): 

cmsun%  cmattach  -b  0.95  -p  16k  -v  256k  -il  wave_3d  & 


or 


cmx'/.  cmattach  -b  0.95  -p  16k  -v  256k  -C  S  wave_3d  & 

3.4  Viewing  the  results 

In  order  to  view  the  results  of  the  3-dimensional  synthetic  seismogram  program,  2-dimensional 
’slices’  through  the  center  of  the  model  may  be: 

1.  output  as  ASCII  files  and  transferred  to  a  VAXstation  for  plotting,  or 

2.  output  to  the  Connection  Machine  frame  buffer  during  runtime. 


3.4.1  Plotting  ASCII  files  on  the  VAXstation 


Program  wave_3d  writes  two  ASCII  files  (one  containing  divergence  and  the  other  containing 
curl)  for  each  2-dimensional  plane  selected,  every  time  step  interval.  The  filenames  are  determined 
by  the  2-dimensional  plane  and  the  time  step: 


output  option 
out-xy 

out.yz 

out_xz 


files 

divxy_N.dat 

curlz_N.dat 

divyz_N.dat 

curlx_N.dat 

divxz_N.dat 

curly_N.dat 


description 

xy  plane  at  z  =  z-dimension/2-1 
yz  plane  at  x  =  x-dimension/2-1 
xz  plane  at  y  =  y-dimension/2-1 


where  N  is  the  time  step. 


The  procedure  for  plotting  ASCII  files  of  divergence  and  curl  on  the  VAXstation  is  as  follows: 


1.  Transfer  the  ASCII  files  from  the  CM  front-end  to  the  RED  VAX  using  ftp: 

RED-3  ftp  134.207.7.12  (refer  to  section  4.4) 

username:  alien 

password:  <<type  password  here>> 

*  get  divxz_200.dat  (to  transfer  divxz-200.dat  from  cmvax  to  RED) 

*  get  curly_200.dat  (to  transfer  curly.200.dat  from  cmvax  to  RED) 

*  ex  (to  exit  ftp) 
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Figure  3:  Example  of  cmattach  command 


Figure  4a  -  Example  of  cinfinRcr  and  emusers  commands  on  the  NRL  cmvax  front-end 
cavax'/.  cmfinger 
Connection  Machine:  BAMBI 


nobody 

CMSUH : 0 

THIRK7E 

Whitcomb 

CMVAX :0 

Connection  Machine: 

GODZILLA 

ancona 

CMSUH : 1 

webb 

THIHK40 

jennings 

CMVAX :1 

cmvax*/,  emusers 

Interface  User 

CM-Idle 

0  vhitcomb 

1  jennings 

:  03 

cmvax/. 

Rot  attached. 

Rot  attached. 
Sequencer  ports  (0) 


Sequencer  ports  (1) 
Rot  attached. 

Rot  attached. 


Status 

ATTACHED  running  “itarlisp" 
ATTACHED  running  "cnattach" 


Figure  4b  -  Example  of  cmfinger  and  emusers  commands  on  the  NRL  cmvax  front-end 
cosun'/,  cmfinger 
Connection  Machine:  BAMBI 


tedwards 

CMSUH : 0 
THIRK76 

aandelbe 

CMVAX :0 

Connection  Machine: 

GODZILLA 

gllen 

CMSUH :1 

uebb 

THIHK40 

nobody 

CMVAX :1 

casun'/,  emusers 

Interface  User 

0  tedwards 

1  alien 

CM-Idle 

-  root 

caenn% 

:47 

Sequencer  ports  (1) 
Rot  attached. 
Sequencer  ports  (0) 


Sequencer  ports  (0  1) 
Rot  attached. 

Rot  attached. 


Status 

ATTACHED  running  "emattach" 
ATTACHED  running  "»aYe_3d" 
DETACHED 


Figure  4c  -  Example  of  cmfinger  and  emusere  commands  on  the  NPAC  cmx  front-end 
car'/,  cmfinger 

Connection  Machine:  THE-CUBE 


nobody  CMX . RPAC . S YR . EDU : 0  Hot  attached. 


Connection  Machine:  S0H-0F-C0BE 


nobody 

root 

nobody 


cnxX  emusers 

Interface  User 
2  root 


CKX.RPAC.SYH.EDU:1  Rot  attached. 

CMX . RPAC . SYR . EDU : 2  Sequencer  ports  (3) 
CMX. IP AC. SYR. EDU: 3  Hot  attached. 

CM1S  Rot  attached. 


CM-Idle  Status 
1:62  ATTACHED  running  "1" 


cmxX 


Figure  4:  Example  of  cmfinger  and  emusers  commands 
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RED-S 


2.  On  the  RED  VAX,  convert  the  ASCII  output  files  to  binary  format  using  the  convert  program: 
RED-8  ©convert 

Enter  CM  (input)  filename:  divxz_200.dat 
Enter  binary  output  filename:  cmxz0200.div 
Enter  scale  factor  (1  for  no  scaling):  1 
FORTRAI  STOP 
RED -8  ©convert 

Enter  CM  (input)  filename:  curly_200.dat 
Enter  binary  output  filename:  cmxz0200.crl 
Enter  scale  factor  (1  for  no  scaling):  1 
FORTRA*  STOP 
RED -8 

3.  Log  onto  the  VAX  Workstation  THISLE  (THISLE  is  on  the  VAX  cluster  so  the  default  data 
directory  RNR2:[FIND.JMA1J  can  be  accessed  from  either  RED  or  THISLE). 

4.  On  THISLE,  run  snap-fix  to  enlarge  the  images,  if  desired  (this  example  will  enlarge  a  64x64 
image  to  256x256): 

THISLE  J$  r  snap -fix 

Enter  input  filename  (fin) :  cmxz0200.div 
Enter  output  filename  (fout):  cbxz0200.div 
Enter  x,z  dimensions  (max  128):  64,64 
Enter  xinx.zinc:  4,4 
FORTRA*  STOP 
THISLE.S 

5.  On  THISLE,  execute  the  command  SNAP  to  send  the  image  to  the  workstation  screen  and/or 
create  an  image  file.  The  example  in  figure  5  creates  an  image  file  named  cbxz0200.img. 

6.  To  get  a  hardcopy  of  the  image  file(s)  created  by  SNAP,  log  back  onto  the  RED  VAX  and  run 
the  command  file  img.  Output  is  to  the  Imagen  laser  printer  (figure  6): 

RED-8  ©img  cbxz0200.img 


3.4.2  Using  the  CM  frame  buffer 

The  CM  graphics  display  system  (GDS)  consists  of  a  frame  buffer  and  a  high-resolution  color 
monitor  (19  inch,  1280  x  1024  pixels),  which  combine  to  display  images  computed  on  the  CM-2. 
The  GDS  can  be  programmed  to  display  any  2-dimensional  grid  of  virtual  processors. 

If  the  variable  oui.fb  is  set  to  1,  2  or  3  in  the  file  xmodel.dat,  images  will  appear  on  the  default 
frame  buffer  at  NRL  or  NPAC  in  ’real-time’  (this  is  obviously  only  useful  if  there  is  someone  at  NRL 
or  NPAC  who  is  watching  the  monitor!). 

At  NRL,  the  four  frame  buffer  locations  are  as  follows: 
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<«  note  that  caps  are  required  for  portions  of  this  program  >>> 


THISL£_$  SNAP 
ENTER  FILEID 
CBXZ 

ISTAT=  0 

ERROR  29  OPENING  FILE  CBXZ. PAR 

ISTAT=  29 

ENTER  MM,NH,ND,DELZ 

64,64,0,50 

ENTER  TIMESTEP 

200 

ENTER  TOP  FILE  EXTENT 
DIV 

EHTER  BOTTOM  FILE  EXTENT  ("H"  FOR  NO) 

CRL 

FILE  OPENED:  CBXZ0200.DIV 
FILE  OPENED:  CBXZ0200.CRL 
ENTER  REDUCE  FACTOR 
1 

1 

ENTER  RANGE  MIN, MAX  (POINTS) 

1.64 

1  64 

ENTER  DEPTH  MIN, MAX  (POINTS) 

1.64 

1  64 

ENTER  CLIPPING  MIN .MAX  VALUES 

-0.01,0.01 

CLIPPING  VALUES =  -9 .9999998E-03  9 . 9999998E-03 

ENTER  ANNOTATION  FOR  TITLE 
«<  hit  return  here  >» 

OUPUT  TO  WORKSTATION  MONITOR?  (Y/N) 

Y  «<  use  the  workstation  mouse  to  manipulate  the  menus  »> 
NEW  COLOR  MAP  ?  (Y/N) 

N 

RESCALE  ?  (Y/N) 

N 

OUTPUT  IMAGE  ?  (Y/N) 

Y 

FORTRAN  STOP 
THISLE_$ 


Figure  5:  Example  of  running  SNAP  on  THISLE 
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Section  along  the  x-z  plane  at  time  step  200  of  3-dimensional 
synthetic  seismogram  calculated  with  program  wave-3d.  The 
model  calculated  is  shown  in  figure  2.  Image  displayed  on  microVAX 
workstation  using  program  SNAPSHOT. 


Figure  6:  Hardcopy  image  on  the  laser  printer 
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CM 

interface 

sequencer 

FB  location 

Bambi 

0 

0 

User  room 

Bambi 

0 

1 

Anderson/Whaley  office 

Godzilla 

1 

0 

Room  with  coffee  maker/fridge 

Godzilla 

1 

1 

Shirron’s  office 

At  NPAC,  the  frame  buffers  are  attached  to  sequencers  0  and  2  of  SON-OF-CUBE. 


4  Guide  to  programming  on  the  Connection  Machine 

The  first  step  in  writing  a  program  to  run  on  the  Connection  Machine  is  to  learn  to  think  in 
parallel  terms.  ’’Learning  to  write  programs  for  parallel  machines  requires  thinking  in  ways  that  are 
quite  different  from  those  demanded  by  sequential  computers.”  (Hillis,  1987) 


4.1  Special  cases  and  boundary  conditions 

The  large  number  of  processors  in  the  Connection  Machine  are  controlled  by  the  front-end 
computer,  which  sends  a  single  instruction  to  all  the  processors  simultaneously.  For  every  instruction, 
each  processor  must  either  execute  or  not  execute.  Special  cases  and  boundary  conditions  are  one 
weakness  of  parallel  processing.  Algorithms  which  contain  special  cases  should  be  redesigned,  if 
possible. 

For  example,  consider  the  calculation  of  boundary  conditions.  On  a  serial  machine,  the  interior 
of  the  model  would  be  calculated  with  one  equation  and  the  edges  of  the  model  with  another.  On  the 
Connection  Machine,  however,  it  is  most  efficient  for  every  processor  to  perform  the  same  equation 
simultaneously.  One  solution  is  to  add  a  damping  term  to  the  wave  equation  to  eliminate  unwanted 
reflections  from  the  edges  of  the  model.  Then  the  same  equation  can  be  used  for  every  processor 
with  only  the  damping  parameter  changing  with  position  on  the  grid  (Charrette,  1987). 


4.2  Timing  considerations 

Another  power  of  programming  on  a  massively  parallel  computer  such  as  the  Connection  Ma¬ 
chine  is  that  increasing  the  size  of  the  model  (within  the  hardware  limits  of  the  computer)  will  not 
change  the  time  required  to  perform  the  parallel  calculations,  although  serial  operations  such  as 
initialization  will  take  longer.  For  example,  wave-3d,  the  3D  synthetic  seismogram  program,  has 
three  main  sections: 


1.  code  to  read  model  parameters  and  perform  parallel  initialization 

2.  code  to  perform  parallel  calculations  within  a  serial  time  loop 

3.  serial  code  for  output  of  results 

It  is  difficult  to  make  an  accurate  comparison  of  runtimes  on  the  CM  and  runtimes  on  the  VAX 
8800.  The  timing  facility  on  the  CM  gives  real  time  and  CM  time  (time  spent  in  parallel  mode) 
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while  the  VAX  timing  facility  gives  real  time  and  CPU  time.  Comparison  of  real  times  is  not  entirely 
valid  because  of  multi-tasking  on  the  computers.  In  this  report,  we  are  comparing  real  time  on  the 
Connection  Machine  with  CPU  time  on  the  VAX;  the  performance  of  the  Connection  Machine  is 
actually  under-estimated. 

A  3D  model  with  dimensions  64x128x64  and  200  time  steps  took  0.5  seconds  (real  time)  for 
initialization  and  1253.6  seconds  (real  time)  total  on  the  32k  Connection  Machine  at  NPAC. 

A  3D  model  with  dimensions  80x80x50  and  200  time  steps  took  400  minutes  (CPU  time)  on  the 
VAX  8800.  The  following  table  summarizes  the  timing  calculations: 

Computer  Number  of  grid  points  seconds  /  time  step 
CM-2  (32k)  524288  (64x128x64)  "6.3  (real  time) 

VAX  8800  320000  (80x80x50)  120  (CPU  time) 

The  Connection  Machine  is  approximately  twenty  times  faster  than  the  VAX  for  a  model  with 
60%  more  points.  Although  models  computed  on  the  Connection  Machine  are  constrained  by  the 
number  of  virtual  processors  available,  it  may  be  possible  to  compute  larger  models  using  the 
DataVault,  a  high  speed  mass  storage  device  attached  to  the  Connection  Machine  (see  section  6.0). 


4.3  Converting  serial  FORTRAN  to  parallel  C* 

This  section  outlines  some  special  features  of  C*  and  some  guidelines  to  assist  a  programmer  in 
developing  and  modifying  C*  code  for  the  Connection  Machine. 

The  original  3-dimensional  synthetic  seismogram  program  was  written  in  FORTRAN  and  devel¬ 
oped  to  run  on  the  VAX,  a  serial  computer.  The  conversion  of  the  serial  FORTRAN  code  to  parallel 
C*  code  was  accomplished  in  four  basic  steps: 

1.  convert  from  FORTRAN  to  standard  C 

2.  determine  portions  of  the  code  which  could  be  run  in  parallel 

3.  convert  from  standard  C  to  parallel  C* 

4.  optimize  the  C*  code 


4.3.1  Convert  from  FORTRAN  to  standard  C 

One  of  the  most  common  coding  errors  in  converting  from  FORTRAN  to  C  is  in  indexing. 
FORTRAN  indexes  from  1  to  N  while  C  indexes  from  0  to  N-l.  For  example,  a  source  located  at  the 
center  of  a  64x64x64  model  would  have  indices  (32,32,32)  in  FORTRAN  while  the  source  located  at 
the  center  of  the  same  model  would  have  indices  (31,31,31)  in  C. 


4.3.2  Determine  portions  of  the  code  which  could  be  run  in  parallel 

The  C*  programming  language  was  designed  specifically  for  the  Connection  Machine.  C*  is  an 
extension  of  the  C  programming  language  that  incorporates  some  of  the  object-oriented  features  of 
C-H-  to  enable  parallel  processing. 
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The  C*  language  encourages  the  programmer  to  think  in  terms  of  groups  of  related  quantities, 
similar  to  the  struct  concept  in  standard  C.  In  C*,  the  keyword  domain  extends  the  struct  concept 
by  allowing  functions,  as  well  as  data  objects,  to  be  associated  members  of  a  group.  Domains  are 
based  on  the  class  concept  from  C++.  Once  a  domain  is  defined  and  member  space  is  allocated, 
each  instance  of  the  domain  can  be  thought  of  as  having  its  own  virtual  processor.  A  program  can 
contain  both  serial  and  parallel  code;  code  is  parallel  only  within  a  member  function  of  a  domain 
or  within  a  selection  statement  which  has  activated  a  domain.  Otherwise,  the  parallel  code  appears 
just  like  serial  code. 

Parallel  variables  are  the  default  when  declared  within  a  member  function  and  within  a  selection 
statement  and  can  be  explicitly  declared  with  the  poly  keyword  within  serial  code.  Serial  variables 
require  the  keyword  mono  if  declared  within  a  parallel  context. 

Many  of  the  computations  in  the  original  synthetic  seismogram  program  involved  large,  3- 
dimensional  FORTRAN  do  loops  where  the  same  computation  was  performed  for  each  node  of 
the  3-dimensional  model.  These  computations  are  relatively  easy  to  think  of  in  a  parallel  sense. 
Instead  of  processing  each  node  of  the  model  in  sequence,  all  nodes  are  processed  in  parallel  by 
placing  the  same  code  on  each  of  the  parallel  processors.  The  first  step  is  to  allocate  a  virtual  pro¬ 
cessor  for  each  grid  point  using  the  nd-grid  library  routine  make-grid  (see  Appendix  A).  A  domain 
with  23  common  variables  and  20  parallel  functions  was  defined  for  the  3D  model  (Figure  7).  Once 
the  parallel  domain  is  defined,  operations  which  require  a  3-dimensional  do  loop  in  FORTRAN  are 
reduced  to  a  single  operation  in  C*  (Figure  8). 


4.3.3  Convert  from  standard  C  to  parallel  C* 


The  following  C*  header  files  are  available  for  parallel  C*  programs: 


<stdio.hs> 

<math.hs> 

<nd.grid.hs> 

<fb.hs> 

<cm/cmtimer.hs> 

<cm/cmfs.hs> 

<cm/cm-file.hs> 


standard  i/o  library  header  file 

math  routines  header  file 

header  for  n-dimensional  NEWS  routines 

header  for  frame  buffer  routines 

CM  timer  functions 

CM  file  structure  (DataVault)  functions 
DataVault  functions 


4.3.4  Optimize  the  C*  code 

Each  processor  on  the  CM-2  has  64k  bits  of  memory  and  this  was  found  to  be  the  major 
constraint  on  the  maximum  size  model  that  could  be  run.  Techniques  used  to  decrease  the  amount 
of  memory  required  by  the  program  and  maximize  model  size  include: 

1.  re-compute  values  rather  than  store  them  in  memory 

2.  break  up  large  functions  into  smaller  functions  in  order  to  reduce  the  amount  of  memory 
required 

There  is  a  ’bug’  in  the  C*  compiler  which  causes  any  floating  point  constant  which  is  not  implicitly 
appended  with  an  / to  be  double  precision.  This  causes  any  calculations  using  that  constant  to  be 
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/*  common  definitions  for  3D  finite  difference  modeling  */ 


#def ine 

XDIM 

64 

#def ine 

YDIM 

64 

#def ine 

ZDIM 

64 

#def ine 

TOTAL 

XDIM* YDIM* ZDIM 

/*  define  the  domain  for  the  parallel  variables  -  27  variables  =  108 
bytes  */ 
domain  model  { 

float  u0,ul,u2; 
float  v0,vl,v2; 
float  v0,srl,w2; 
float  clamb; 
float  emu; 
float  clp2mu; 
float  rho; 

float  txxdx.txydy.txzdz; 
float  txydx.tyydy.tyzdz; 
float  txzdx,tyzdy,tzzdz; 
float  alpha; 

float  div , curlx , curly , curlz ; 

float  alpha_fct(int  indx); 

void  boundary (void) ; 

void  calc_adisp(float  *adisp) ; 

void  calc_bdisp(float  *bdisp); 

void  calc_cdisp(float  *cdisp) ; 

void  calc_displ(void) ; 

void  calc_txxdx(void) ; 

void  calc_txydy(void) ; 

void  calc.txzdz(void) ; 

void  calc_txydx(void) ; 

void  calc_tyydy(void) ; 

void  calc_tyzdz(void) ; 

void  calc_txzdx(void) ; 

void  calc_tyzdy(void) ; 

void  calc_tzzdz(void) ; 

void  divcurl(void) ; 

void  initialize (void) ; 

void  output.fb(void) ; 

void  stress (void) ; 

void  update(void) ; 

>  point [TOTAL] ; 


Figure  7:  Declaration  of  domain  for  3D  synthetic  seismograms 
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FORTRAN  subroutine  update  -  update  arrays  after  completion 
of  a  time  step 

subroutine  update 

include  ’dspall.blk’ 

common/values/  ltmout , irout , iy out , izout , nr , ny , nz , nt , nt 1 

do  300  iz  =  1,  nz 

do  200  iy  =  1 ,  ny 

do  100  ix  =  1,  nr 

uold(iz,iy,iz)  =  u(ix,iy,iz) 
vold(ix,iy ,iz)  =  v(ix,iy,iz) 
wold(ix,iy , iz)  =  w(ix,iy,iz) 

u(ix,iy,iz)  =  unew(ix,iy,iz) 
v(ix,iy,iz)  =  vnew(ix,iy ,iz) 
w(ix,iy,iz)  =  wnew(ix,iy,iz) 
continue 
continue 
cont inue 
return 
end 


/*  */ 

/*  C*  function  update  -  update  arrays  after  completion  of  */ 

/*  a  time  step  */ 

/*  */ 

#include  "wave_3d.hs" 


void  model: :update(void)  /*  function  type  is  domain  model  */ 

{ 

uold  =  u; 
void  =  v; 
wold  =  w; 

u  =  unew; 
v  =  vnew; 
h  =  unew; 

> 


c======= 

c 

c 

c===- -== 

c 

c 


c 


100 

200 

300 


Figure  8:  Example  of  serial  FORTRAN  vs  parallel  C* 
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done  in  double  precision,  thus  wasting  valuable  memory  resources.  To  avoid  this  problem,  append 
an  /  to  all  floating  point  constants  in  the  C*  program,  for  example: 

♦define  PI  3.141591  (instead  of  ♦define  PI  3.14159) 

alpha  =  sin( ( (ix-(nx-ITEL) )/(ITEL/2 . Of ) )*PI/2. Of )*(0 . 11 ) ; 

To  detect  lines  containing  values  which  would  cause  double  precision  computations,  use  the  -w 
flag  with  the  compile  command: 

cm x'l  cs  -O  -w  -o  progname  progname.es  -lnrlcstar  & 


There  is  a  timing  facility  in  C*  for  accumulating  and  reporting  elapsed  real  time  and  CM  time. 
To  use  the  timing  facility,  include  the  header  file: 

#includ«  <cm/cmtimer .hs> 


The  timer  functions  used  in  the  current  program  are: 

CM-start_timer(l)  begin  accumulating  timing  information 
CM_stop_timer(l)  stop  accumulating  and  print  timing  information 

CM_reset_timer()  reset  the  clock 

These  functions  were  used  to  time  the  initialization  stage  of  the  program  and  the  total  run-time. 
T1  e  timer  functions  can  be  implemented  as  C  preprocessor  macros  as  follows: 

♦define  TIMER -01  CM-start-timsr(l) 

♦define  TIMER-RESET  CM-reset-tiaerO 


♦define  TIMER-OFF  CH-stop.timer(l) 


4.4  File  transfer  via  ftp 

The  internet  connection  to  NRL  is  very  slow.  For  this  reason,  most  of  the  program  coding  and 
editing  was  done  on  the  Red  VAX  at  WHOI.  Programs  were  then  transferred  to  the  CM  front-end 
over  the  network  using  file  transfer  protocol  (FTP),  illustrated  in  Figure  9. 

1.  ftp  to  connect  to  the  CM  front-end 

RED-t  ftp  134.207.7.12  <<<  to  connect  to  the  NRL  VAX  >>> 

RED-1  ftp  134.207.7.4  <<<  to  connect  to  the  NRL  SUN  >>> 

RED .1  ftp  cmx.npac.syr.edu  <<<  to  connect  to  the  NPAC  VAX  >>> 
cmx’/i  ftp  red.whoi.edu  <<<  to  connect  to  the  WHOI  RED  VAX  >>> 

2.  cd  to  the  desired  directory  if  files  are  in  a  directory  other  than  HOME 

3.  use  get  to  copy  files  FROM  the  CM  front-end,  use  put  to  copy  files  TO  the  CM  front-end 
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4.  ex  to  exit  ftp  from  the  WHOI  VAX,  (bye  to  exit  from  the  UNIX  front-ends) 


4.5  Debugging  C*  Code 

Debugging  the  parallel  portions  of  a  C*  program  can  be  difficult.  It  is  generally  best  to  start 
with  a  simple  model  and  print  out  intermediate  values  for  selected  processors.  If  a  printf  statement 
is  used  within  a  parallel  context,  it  must  be  immediately  followed  with: 

iilush(stdout) ; 

This  will  ensure  that  all  buffered  data  is  sent  to  the  output  stream  before  parallel  commands  are 
executed. 

Section  5.3.3  of  the  Connection  Machine  System  C*  User’s  Guide  explains  how  to  run  C*  pro¬ 
grams  with  the  dbx  debugger. 


4.6  Compile  and  link  CM  programs 

Source  code  to  be  run  on  the  CM  should  have  the  file  extent  .cs  and  header  files  should  have 
the  extent  .hs.  Figure  10  illustrates  compilation  of  wave_3d  on  the  NRL  and  NPAC  front-end 
computers. 

To  compile  and  link  a  C*  program  to  run  on  the  NRL  Connection  Machine  (the  NRL  C*  library 
is  installed  as  nrlcstar): 

cmvax'/.  cs  -O  -w  -o  progname  prognaxne.es  -lnrlcstar  tc 


where: 

-O 

-w 

-o  progname 

instead  of  a.out 

progname.es 

-lnrlcstar 

-lemfb 

k 


sets  the  flag  for  the  optimizer 

causes  double  precision  calculations  to  be  flagged 

causes  the  executable  file  to  be  named  progname 

is  the  C*  source  code  file 
links  the  program  with  the  NRL  C*  library 
links  the  program  with  the  frame  buffer  library 
the  UNIX  command  for  running  in  the  background 


To  compile  and  link  a  C*  program  to  run  on  the  NPAC  Connection  Machine  (the  NRL  C*  library 
routines  are  local): 

cmx'/.  cs  -O  -w  -o  progname  progname.es  nd-grid.o  fb.o  & 


4.7  Running  a  C*  program  on  the  Connection  Machine 

To  run  a  program  on  the  NRL  Connection  Machine: 
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itting  n  file  from  NPAC  imx  using  file  transfer  protocol  (ftp) 
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Figure  9:  Example  of  using  file  transfer  protocol  (ftp) 


Figure  10a  -  Example  of  compiling  wave_3d  on  the  NRL  cmsun  front-end 
cmsun*/.  cs  -O  -w  -o  wave_3d  waveJJd.es  -lnrlcstar  & 

[1]  18467 

wavo_3d.es:  Co-compiling  =>  wave_3d..c 
Co  Compiler  Version  5.0.21 

"wave_3d.es",  line  98:  warning:  Underflow  from  floating-point  literal  0.0. 
"wave_3d.es",  line  98:  warning:  Underflow  from  floating-point  literal  0.0. 
"wave_3d.es",  line  100:  warning:  Underflow  from  floating-point  literal  0.0 
"wave_3d.es",  line  977:  warning:  Underflow  from  floating-point  literal  0.0 
wave_3d.es:  C-compiling  wave_3d..c 
Linking 

[1]  Done  cs  -0  -w  -o  wave_3d  wave_3d.es  -lnrlcstar 

cmsun*/. 


Figure  10b  -  Example  of  compiling  wave_3d  on  the  NPAC  cmx  front-end 
cmx*/.  cs  -O  -w  -o  waveJJd  wave-3d.es  nd-grid.o  fb.o  & 

Cl]  5996 

cmx*/.  wave_3d.es:  C*- compiling  =>  wave_3d..c 

C*  Compiler  Version  5.0 

wave_3d . cs :  C-compiling  wave_3d . . c 

Linking 

Cl]  Done  cs  -0  -w  -o  wave_3d  wave_3d.es  nd_grid.o  fb.o 

cmx*/. 


Figure  10:  Compilation  of  waveJJd 


cmeun%  cmattach  -b  0.95  -p  nprocs  -v  nvprocs  -i  interface  progna me  >  output  Jile  & 


To  run  a  program  on  the  NPAC  Connection  Machine: 

cmx*/,  cmattach  -b  0.95  -v  nvprocs  -p  nprocs  -C  CM -name  progname  >  output  Tile  & 
where: 

causes  95%  of  the  memory  per  virtual  processor  (vp)  to  be  available  (default 
is  0.75) 

attach  nprocs  CM  processors;  legal  values  for  nprocs:  4k  or  8k  for  BAMBI, 
(default  is  4k)  8k  or  16k  for  GODZILLA  (default  8k)  8k,  16k  or  32k  for 
SON-OF-CUBE  (default  8k) 

configure  the  CM  to  have  nvprocs  virtual  processors  (must  be  at  least  the 
number  of  physical  processors  allocated) 

attaches  the  CM  via  interface  0  (BAMBI)  or  interface  1  (GODZILLA)  on 
NRL  CM 

attaches  to  the  CM  named  ’CM_name’  (use  -C  S  for  SON-OF-CUBE  at 
NPAC) 

your  executable  C*  program 
file  containing  log  of  batch  session 
causes  program  to  run  in  the  background 

The  Connection  Machine  System  C*  User’s  Guide  (section  4.3)  includes  instructions  for  executing 
C*  programs  interactively. 

4.8  Re-configuring  wave_3d 

The  wave-3d  program  is  currently  configured  for  models  with  dimensions  64x64x64,  which  is 
the  largest  model  that  will  run  on  a  16k  Connection  Machine.  The  32k  Connection  Machine  at 
NPAC  can  accomodate  a  model  with  dimensions  64x128x64.  To  change  the  model  dimensions  in 
the  wave -3d  program: 

1.  edit  the  header  file  wave-3d.hs  and  modify  any  or  all  of  the  following  lines  (note  that  model 
dimensions  MUST  be  a  power  of  2!): 

•define  XDIM  64 

•define  YDIM  64 

•define  ZDXM  64 

2.  re-compile  the  wave_3d  program  (refer  to  section  4.6) 
cox'/,  cs  -O  -w  -o  wave_3d  wave-3d.cs  nd-grid.o  fb.o  Ic 

3.  edit  xmodel.dat  to  incorporate  the  new  model  parameters 


-b  0.95 
-p  nprocs 

-v  nvprocs 

•i  interface 

-C  CM  .name 

progname 
output  .file 
& 
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4.  attach  and  run  the  program,  making  the  appropriate  changes  to  the  cmattach  parameters 
(refer  to  section  4.7) 

cmx'/.  cmattach  -b  0.95  -p  32k  -v  512k  -C  S  wave_3d  >  wave_3d.out  & 


5  Frame  Buffer 


The  Connection  Machine  graphic  display  system  consists  of  a  framebuffer  and  a  high  resolution 
color  monitor.  The  framebuffer  contains  a  display  memory  capable  of  storing  an  image  2048  by 
1024  pixels  with  three  eight-bit  color  buffers  and  a  four-bit  overlay  buffer.  The  color  monitor  has  a 
resolution  of  1280  by  1024  pixels.  Since  the  display  memory  can  store  more  pixels  than  the  monitor 
can  display,  pan  and  zoom  operations  are  supported. 

The  CM  graphic  display  system  allows  each  virtual  processor  to  be  mapped  to  a  pixel  on  the  color 
monitor.  As  values  from  the  virtual  processors  are  written  to  the  display  memory,  the  image  on  the 
monitor  changes.  For  the  3D  synthetic  seismograms,  a  2-dimensional  plane  is  selected  (for  example, 
the  xz  plane)  and  ’slices’  through  the  model  are  displayed  at  selected  time  steps  in  ’real-time’. 

The  framebuffer  is  normally  accessed  via  C/PARIS  (Connection  Machine  PARallel  Instruction 
Set)  calls.  The  wave -3d  program  uses  the  NRL  C*  library  fb  which  is  a  collection  of  functions 
used  to  display  pixels  on  the  framebuffer  (see  Appendix  A).  To  use  these  functions  in  a  C*  program, 
include  the  header  file: 

# include  "tb.cn" 

and  link  with  the  library:  -lnrlcstar  (if  using  NRL  Connection  Machine)  or  fb.o  (if  using  NPAC 
Connection  Machine). 


6  DataVault 


The  DataVault  is  a  5-Gigabyte  mass  storage  system  which  provides  a  data  parallel  file  system 
for  the  Connection  Machine.  Data  can  be  transferred  to  and  from  the  CM-2  in  parallel  at  high  speed 
(40  megabytes  per  second)  or  it  can  be  read  and  written  directly  from  the  front-end  without  using 
the  CM-2. 

The  amount  of  memory  per  processor  and  the  number  of  physical  processors  on  the  Connection 
Machine  constrains  the  maximum  size  of  the  3D  model  that  can  be  computed. 

In  order  to  calculate  3D  synthetic  seismograms  for  models  with  dimensions  of  128x128x128 
(actual  dimensions  124x128x124  with  overlap),  the  model  can  be  divided  into  four  segments  with  an 
overlap  of  4  grid  points  (figure  11).  Each  segment  can  be  copied  to  a  temporary  file  on  the  DataVault 
at  each  time  step.  Preliminary  calculations  indicate  that  approximately  2  seconds  per  time  step  are 
required  to  write  the  23  parallel  domain  variables  to  four  temporary  files  on  the  DataVault. 
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Three-dimensional  model  divided  into  four  segments;  each  segment  has  di¬ 
mensions  64x128x64.  The  segments  overlap  four  grid  points  in  the  x-  and  z- 
directions,  giving  total  model  dimensions  of  124x128x124.  At  every  time  step, 
the  variables  calculated  for  each  segment  will  be  stored  on  the  DataVault. 


Figure  1 1 :  Segmentation  of  large  model  using  DataVault 
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Appendices 


A  NRL  cstar  library 


The  NRL  cstar  library  consists  of  functions  written  to  provide  a  simple  interface  to  many  of  the 
Connection  Machine  features  that  can  only  be  accessed  directly  via  C/Paris.  The  library  modules 
were  written  to  solve  specific  problems  encountered  by  users  at  NRL.  A  full  description  of  the  library 
can  be  found  in  the  README  file  located  in  the  nrl-cstar-lib  directory  on  the  NRL  and  NPAC 
front-end  computers  (cmvax,  cmsun  and  cmx). 

The  NRL  C*  library  resides  on  the  NRL  front-end  computers  in  the  directory: 

/ca-lights/library /nrl-cstar-lib 

On  the  NPAC  front-end  computer,  the  NRL  C*  library  is  located  in  the  directory: 
/public/cm-lights/nrl-cstar-lib 


A.l  N-dimensional  grid 

The  finite  differences  method  uses  many  ’nearest-neighbor’  calculations.  Although  C*  supports 
a  NEWS  package  to  perform  nearest  neighbor  calculations  in  two  dimensions,  this  package  currently 
does  not  support  three-dimensional  NEWS  communications.  The  nd-grid  module  of  the  NRL  C* 
library  provides  fast  nearest  neighbor  communications  on  N-dimensional  grids  (up  to  8  dimensions). 

For  more  information  on  the  NRL  C*  n-dimensional  grid  library(figure  12): 

cuxV,  cd  /public/cm-lights/nrl-cstar-lib/nd-grid 
cmx'/,  more  README 


A. 2  Framebuffer 

The  NRL  C*  library  module  fb  is  a  simple  interface  to  the  high  speed,  high  resolution  CM 
graphics  display  system.  This  module  consists  of  functions  to  initialize  the  framebuffer,  set  up  a 
color  table,  and  display  pixels  on  the  color  monitor. 

For  more  information  on  the  NRL  C*  framebuffer  library  (figure  13): 

cbx%  cd  /public/cm-lights/nrl-cstar-lib/fb 
cmx%  more  README 


B  cmusers  group 


NRL  supports  a  users’  forum  for  exchanging  questions,  ideas  and  experiences  among  program¬ 
mers  using  the  Connection  Machine  Facility.  Mail  can  be  directed  to  cmusers  and  questions  usually 
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Author:  Robort  Wholoy 

E— no  1 1 :  wholoy*nr l-c«( .arpo 

US  Mil:  Robort  Wholoy 

Naval  Roooareh  Lab 
Coda  5136 
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Author : 
E-mo i I : 
US  moil: 


Phone ; 

Af  f i I iote: 


Robert  Who  ley 
whaleyOnr l-cmf .orpo 
Robert  Who  ley 
Naval  Reseorch  Lob 
Code  5150 

Washington,  D.C.  20375-5000 

(202)  404-7019 

TMC  (Site  Rep  at  NRL) 


Purpose : 

These  functions  are  used  to  display  pixels  on  the  framebuffer. 

Usage : 

The  framebuffer  is  first  initialized,  then  the  color  map  can  be  set, 
pixels  displayed  and  finally  the  framebuffer  released. 

Syniox: 

void  ini t_f rome_buf f er( i nt  x_size,  Int  y_size); 

Examp  I e : 

ini t_f rame_buf fer(512,  512); 

Attach  the  fromebuffer  and  configure  it  so  that  there  are  512  pixels 
horizontally  and  512  pixels  vertically. 


Syntax : 

set_color( int  color_ld,  Int  red,  Int  green,  int  blue); 

Example: 

set_color(1O0,  128.  0,  128); 

Set  the  color  map  so  that  the  color  100  Is  dork  purple.  The  default 
color  mop  is  gray  scale  from  0  (black),  to  255  (white).  255  colors  can 
be  set.  Red,  green,  ond  blue  con  eoch  range  from  0  to  255  giving  a  total 
of  16  million  possible  shades.  Color  0  is  the  background  color,  so  it  is 
usually  a  good  idea  to  leave  it  block  (red-0,  green-0,  blue-0). 

'  Example: 

for  ( i— 0;  i <25;  i++)  set.col or( i+200.  1*10,  0.  250  -  1*10); 

Set  the  color  map  values  from  200  through  224  to  range  from  bright 
blue  (value  200)  to  bright  red  (value  224)  with  progressive  shades 
of  purple  between. 


Syntax; 

plot_f rom_grld(uneigned  char  eolor); 

Examp I e : 

unsigned  char  heat: 
plot_f ronugrid(heat) ; 

Each  selected  processor  updates  a  single  pixel  to  the  color  specified  by 
the  variable  *heat*.  The  pixels  updated  are  the  ones  with  the  same  x,  y 
coordinates  as  the  processor.  It  Is  an  error  to  call  this  routine  if 
the  processors  are  not  arranged  as  a  2— d  grid. 


Syntax: 

P I ot_x«jr (unsigned  short  x,  unsigned  short  y,  unsigned  char  color); 


unsigned  short  plx_x,  plx_y; 
unsigned  char  ptx_color; 
plot.x^yCP**-**  p!x«y.  plx-color); 


Eoch  selected  processor  sets  the  pixel  at  'plx_x',  'plx_y*  to  the 
color  'plx^color* . 


Syntax: 


void  release_frame_buffer(vold); 


Example: 

r e I eaee_f r ome_bu f f e r ( ) ; 


Frees  the  framebuffer  to  be  used  by  other 


Figure  13:  README  file  for  fb  library 
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receive  quick  and  thorough  responses.  Since  all  members  of  the  cmusers  group  receive  the  mail, 
you  can  learn  about  problems/bugs  with  new  compilers,  upcoming  classes  and  talks,  innovative 
algorithms,  etc.  To  join  this  forum,  simply  mail  your  request  to  be  added  to  the  mailing  list  to 
cmusers. 

To  get  help  on  the  mail  facility: 

CBsuny.  man  mail 

NPAC  supports  an  online  consultant  service.  To  get  help,  send  mail  to  consult. 
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C  Structure  diagram  for  wave_3d.cs 


setup.fb 


input_defs 


wave_3d 
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D  Source  Code 


D.l  wave_3d  program 


/*  program  save_3d  */ 

/*  */ 

/ *  Programmer:  Julie  Allen  */ 

/*  Date:  28  February  1989  */ 

/*  +/ 

/*  */ 

/*  Copyright  (c)  1989  */ 

/•  Information  Processing  and  Communications  Laboratory  */ 

/*  Woods  Hole  Oceanographic  Institution  */ 

/*  All  rights  reserved.  */ 

/+  */ 

/*  This  material  cannot  be  distributed  or  sold  without  */ 

/*  prior  permission  of  the  author(s).  */ 

/*  */ 

/ *  sssrsss:ss:::::s;=ssssssssscsssczsss3:ss3ss:ssss:s:rsscssss  */ 

/*  */ 

/*  */.cs  -0  ~o  save_3d  save_3d.cs  nd_grid.o  */ 

A  '/.cmattach  -b  0.95  -p  16k  -v  256k  -C  S  save_3d  k  */ 

A  */ 

A  */ 

•include  <stdio . hs> 

t include  <cm/cmtimer.hs> 

•include  <math . hs> 

•include  <str ing . h> 

•include  <ctype . h> 

•include  "nd_grid.hs" 

•include  Mfb.hs" 

•include  "save_3d.hs" 


void  mainO 
int  i; 

A  function  declarations  */ 
void  input.defs(void) ; 
void  hetiso(void) ; 

void  make_grid(int  dO,  int  dl,  int  d2); 
void  setup.fb(void) ; 

printf(”\n  **********  Program  save_3d  **********\n\n") ; 
printf ("3D  finite  difference  modeling  programW); 
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print! (“Running  on  Thinking  Machines  CM/2  Parallel  Processor\n") ; 
print! ("Array  size:  */.d  XDIM:  */.d  YDIM:  */.d  ZDIM:  V,d\n" , 

TOTAL,  XDIM,  YDIM,  ZDIM) ; 

/*  define  the  3D  grid  */ 

[domain  model], {  make_grid(XDIM, YDIM, ZDIM) ;> 

/*  setup  the  Irame  bui!er  and  color  table  */ 
setup_!b() ; 

/*  input  model  deiinition  and  parameters  */ 

TIMER.OB;  /*  start  timer  */ 

(void)  input_de!s() ; 

printlO'Time  lor  function  input _de!s :\n") ; 

TIMER_OFF;  /*  stop  timer  */ 

TIMER. 01;  /*  start  timer  */ 

/*  set  up  initialize  conditions  */ 

[domain  model] . {initialize () ;> 
printlO'Time  lor  initialization:^"); 

TIMER. OFF;  /*  stop  timer  */ 

TIMER.OI;  /*  start  timer  */ 

/*  perform  hetergeneous/isotropic  finite  difference  calculations  */ 
(void)  hetisoO; 


printlO'Time  1 or  rest  of  program :\n") ; 

TIMER.OFF ;  /*  stop  timer  */ 

printfO'End  ol  program  wave_3d. cs\n") ; 

> 

/*  */ 

/*  Function  alpha.lct  */ 

/*  */ 

float  model: : alpha.lct (int  indz) 

if (indz  <  ITEL) 

retura(sin(((ITEL-indx-l)/(ITEL/2))*PI/2.f )*(0. iOf )) ; 

else 

retum(sin( ((indx-(ZDIH-ITEL))/(ITEL/2))*PI/2.f)*(0. IOf )) ; 


/* 

/* 

/* 


Function  boundary 


====  */ 

*/ 
*/ 
*/ 
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void  model: : boundary (void) 

{ 

int  ix,iy,iz; 

Iloat  alphx  =  0.01,  alphy  =  0.01,  alphz  =  0.01; 

alpha  =  0.01; 

iz  =  this_coordinate(0) ; 
iy  =  this_coordinate(l) ; 
iz  =  this_coordinate(2) ; 

il(iz  <  ITEL  II  iz  >=  (ZDIM  -  RTEL) ) 

il (iz  <  ITEL/2  II  iz  >=  (ZDIM  -  ITEL/2)) 
alphz  =  0 . 101 ; 

else 

alphz  ■  alpha.lct(iz) ; 

il (iy  <  ITEL  II  iy  >=  (YDIM  -  ITEL)) 

il (iy  <  ITEL/2  II  iy  >=  (YDIM  -  ITEL/2)) 
alphy  *  0.101; 

else 

alphy  =  alpha_lct(iy) ; 

il(ix  <  ITEL  II  iz  >*  (XDIM  -  ITEL)) 

il(ix  <  ITEL/2  II  iz  >*  (XDIM  -  ITEL/2)) 
alphz  *  0.101 ; 

else 

alphz  *  alpha_lct ( iz) ; 

alpha  =  (alphz  >  alphy)  ?  alphz  :  alphy; 
alpha  =  (alphz  >  alpha)  ?  alphz  :  alpha; 

> 

/* 

/*  Function  calc.adisp 

/* 

void  model: :calc_adisp(lloat  *adi*p) 

< 

Iloat  tn0u,tnlu,ta2u,tn0v,tn0*; 

Iloat  tpOu,tpOv,tplv,tpOw; 

/*  set  up  temporary  'next'  values  */ 

tnOu  *  next(O.fcul); 

tnlu  =  nezt(l.Iul); 

tn2u  *  next(2,*ul); 

tnOv  *  next(O.Jtvl) ; 


*=  */ 

*/ 

*/ 

*/ 

=  ♦/ 
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tnOw  =  next(0,fcwl); 


/*  set  up  temporary  ’prev’  values  */ 

tpOu  =  prev(O.Jtul) ; 

tpOv  =  prev(0,*vl); 

tplv  =  prev(l,fcvl); 

tpOu  =  prev(0,*ul); 

/*  term  1:  differences  placed  in  'a'  array 

consists  of  difference  terms  used  in  the  d/dx  stress  term  */ 

*(adisp+0)  =  (tnOu  -  ul)  *  xlamb; 

*(adisp+l)  =  (ul  -  tpOu)  *  xlamb; 

*(adisp+2)  =  (vl  -  tplv)  *  xylamb; 

*(adisp+3)  =  (tpOv  -  prev(l ,*tpOv) )  *  xylamb; 

*(adisp+4)  -  (el  -  prev(2,*vl))  *  xzlamb; 

*(adisp+5)  =  (tpOu  -  prev(2 ,*tpOw) )  *  xzlamb; 

*(adisp+6)  =  (next (2,*tnOu)  -  tnOu)  *  xzlamb; 

*(adisp+7)  ■  (tn2u  -  ul)  *  xzlamb; 

*(adisp+8)  =  (tnOu  -  ul)  *  xlamb; 

*(adisp+9)  =  (ul  -  tpOu)  *  xlamb; 

*(adisp+10)  =  (next(l .fctnOu)  -  tnOu)  *  xylamb; 

*(adisp+ll)  *  (tnlu  -  ul)  *  xylamb; 

*(adisp+12)  *  (tnOv  -  vl)  *  xlamb; 

*(adisp+13)  *  (vl  -  tpOv)  *  xlamb; 

> 

/*  EZSIlSSSSSISSKHSERElZESKaEUEESSKISSSSSCCHnCSKCCCECEISe  */ 

/*  */ 

/*  Function  calc.bdisp  */ 

/*  */ 

void  model :: calc.bdisp(float  *bdisp) 

float  tnOu, tnlu, tnlv,tn2v, tnlu; 
float  tplu,tp0v,tplv,tplu,tp2v; 

/*  set  up  temporary  'next*  values  */ 

tnOu  =  next(O.kul) ; 

tnlu  *  next(l,tul); 

tnlv  =  next(l,fcvl); 

tn2v  =  next(2,fcvl); 

tnlu  =  next(l.kul); 

/*  set  up  temporary  ’prev*  values  */ 
tplu  =  prev(l,*ul); 
tpOv  =  prev(0,*vl); 
tplv  =  prev(l.Avl); 
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tpls  =  prev(l.ksl); 
tp2s  =  prev(2,Zsl); 

/*  term  2:  differences  placed  in  ’b’  array 

terms  used  in  the  d/dy  stress  term  */ 
*(bdisp+0)  =  (next(l ,*tnOu)  -  tnlu)  *  rylamb; 
*(bdisp+l)  =  (tnOu  -  ul)  *  xylamb; 

*(bdisp+2)  =  (tnlv  -  vl)  *  ylamb; 

*(bdisp+3)  =  (vl  -  tplv)  *  ylamb; 

*(bdisp+4)  =  (tnlu  -  prev(2,*tnls))  *  yzlamb; 
*(bdisp+5)  =  (si  -  tp2s)  *  yzlamb; 

*(bdisp+6)  =  (tnlu  -  ul)  *  ylamb; 

*(bdisp+7)  =  (ul  -  tplu)  *  ylamb; 

*(bdisp+8)  =  (vl  -  tpOv)  *  xylamb; 

*(bdisp+9)  =  (tplv  -  prev(l .fctpOv))  *  xylamb; 
*(bdisp+10)  =  (tn2v  -  vl)  *  yzlamb; 

*(bdisp+ll)  =  (next(2,*tplv)  -  tplv)  *  yzlamb; 
*(bdisp+12)  =  (tnls  -  si)  *  ylamb; 

*(bdisp+13)  =  (si  -  tpls)  *  ylamb; 


/*  SSSSSSSSSI*SSS*SS*SSSSSSSS*»*SS**S»»««*XM****IS»IH»»S*S1****  */ 
/*  */ 

/*  Function  calc.cdisp  */ 

/*  */ 

/*  sssssssssssssssssssssssscsrssssssssssBSsssssssstssssssssssssxs  */ 


void  model; : calc_cdisp(f loat  *cdisp) 

float  tn0u.tn2u.tn2v, tnls ,tn2s; 
float  tp2u,tplv,tp2v,tp0s,tp2s; 

/*  set  up  temporary  ’next'  values  */ 

tnOu  =  next(O.hul); 

tn2u  =  next(2,tul); 

tn2v  *  next(2,*vl); 

tnls  -  next(l.ksl); 

tn2s  =  next(2,*sl); 

/*  set  up  temporary  'prev*  values  «■/ 

tp2u  *  prev(2,*ul); 

tplv  *  prev(l,*vl); 

tp2v  =  prev(2,*vl); 

tpOs  *  prev(O.hsl); 

tp2s  *  prev(2,*sl); 

/*  term  3:  differences  placed  in  'c'  array 

consists  of  difference  terms  used  in  the  d/d z  stress  term  */ 
*(cdisp+0)  *  (next(2,ttn0u)  -  tn2u)  *  xzlamb; 

*(cdisp+l)  «  (tnOu  -  ul)  *  xzlamb; 
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*(cdiap+2)  s  (tn2v  -  naxt(2,*tplv))  *  yzlaab; 
*(cdisp+3)  =  (vl  -  tplv)  *  yzlaab; 

*(cdisp+4)  *  (tn2a  -  ul)  *  zlaab; 

*(cdisp+S)  *  (ul  -  tp2u)  *  zlaab; 

*(cdisp+6)  =  (tn2u  -  ul)  *  zlaab; 

*(cdiap+7)  *  (ul  -  tp2u)  *  zlaab; 

*(cdiap+8)  =  (ul  -  tpOu)  *  xzlaab; 

♦(cdiap+9)  *  (tp2u  -  prav(2,*tp0a))  *  xzlaab; 
*(cdisp+10)  =  (tn2v  -  vl)  *  zlaab; 

*(cdiap+ll)  =  (vl  -  tp2v)  *  zlaab; 

*(cdisp+12)  =  (tula  -  al)  *  yzlaab; 

*(cdisp+13)  =  (prav(2,*tnlu)  -  tp2a)  *  yzlaab; 

> 


A  ♦/ 

A  Function  calc_displ  */ 

A  a/ 

void  nodal: :calc_displ( void) 

{ 

int  ix.iy.iz; 

float  alphl , alph2 , alph3 ; 


ix  *  thi»_coordinate(0) ; 
iy  =  this_coordinata(l) ; 
iz  =  thia_coordinata(2) ; 

if  ((ix  <1  II  ix  >«=  XDIM)  II  (iy  <  1  II  iy  >*  YDIH)  I  I 
(iz  <  1  II  iz  >*  ZD IN)) 

» 

alsa 

alphl  =  2.f  -  alpha  *  alpha; 
alph2  =  l.f  -  2.f  a  alpha; 

alph3  =  l.f  +  2.f  a  alpha; 

uO  *  (alphl  *  ul  -  alph2  *  u2  +  txxdx  +  txydy  +  txzdz)/alph3; 

vO  «  (alphl  *  vl  -  alph2  •  v2  +  txydx  ♦  tyydy  +  tyzdz)/alph3; 

aO  *  (alphl  *  al  -  alph2  a  *2  +  txzdx  +  tyzdy  +  tzzdz)/alph3; 

> 


> 

/a  n<iitsMMn»»i»i>iniiini>iniiiii>mninKe:H»i»it  */ 

A  */ 

/a  Function  calc.txxdx  */ 

/a  a/ 

/*  sis>»»:KinniiM»HBK«Knii>iKi»n«cHi>:i»>ttnsis  a/ 
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void  model: :calc_txxdx(void) 

{ 

int  ix,iy,iz; 
float  alp2ul , alp2u2 ; 

float  amuxz 1 , amuxz2 , amuxz3 , amuxz 4 , amuxzS ,  amuxz 6 ; 

float  amuxy 1 , aauxy2 , amuxy 3 , amuxy4 , amuxy 5 , amuxy 6 ; 

float  amuyzl , amuyz2 , amuyz3 , amuyz4, amuyzS , amuyzS ; 

float  amuyz7,amuyz8; 

float  abigul ,abigu2; 

float  alambl , alamb2 ; 

float  temp; 

float  tn0,tnl,tn2,tp0,tpl,tp2; 
float  adisp[14]; 

/*  first  get  the  displacement  difference  terms  */ 
calc.adispCadisp) ; 

/*  set  up  temporary  values  */ 
tnO  =  next(0,fcclp2mu) ; 
tpO  =  prev(0,fcclp2mu) ; 

/*  find  the  lambda* 2mu  values  for  the  txx.tyy.tzz  positions  */ 
alp2ul  «  (clp2mu  +  tn0)/2.f ; 
alp2u2  *  (clp2mu  *  tpO)/2.i; 

/*  set  up  temporary  values  */ 

tnO  a  next (0 ,tcmu) ; 

tnl  a  next (1 .Jfccmu) ; 

tn2  a  next(2,ftcmu) ; 

tpO  a  prev(O.Jtcmu) ; 

tpl  a  prev(l ,*cmu) ; 

tp2  a  prev(2,Jtcmu)  ; 

/*  next,  the  mu  value  for  the  txz  positions  */ 
amuxz 1  a  (emu  *  tn2)/2.f; 

amuxz2  =  (next(0,kcmu)  *  next(2,*tn0))/2.f ; 
amuxz3  a  (emu  *  prev(2,temu))/2.f ; 
amuxz4  *  (next(0,tcmu)  +  prev(2,*tn0))/2.f ; 
amuxzS  *  (prev(0,*emu)  +  next(2,*tp0))/2.f ; 
amuxzfi  «  (prev(0, Jtcmu)  ♦  prev(2,*tp0))/2.f; 

/•  non,  the  mu  for  the  txy  positions  */ 

amuxy 1  *  (emu  *  next(l,*cmu))/2.f; 

amuxy 2  a  (next(0,tcmu)  *  next(l ,fctn0))/2.f ; 

amuxy3  *  (emu  +  prev(l,*cmu))/2.f ; 

amuxy 4  *  (next(O.tcmu)  +  prev(l,fctn0))/2.f ; 

amuxy 5  »  (prev(0 , Jtcmu)  +  next (1  ,*tp0) )/2 . f ; 

amuxyfi  «  (prev(0,*cmu)  prev(l ,*tp0))/2.f ; 
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/*  finally,  the  nu  for  the  tyz  positions  */ 


temp  =  next(l ,fctn0) ; 

amuyzl  =  (emu  +  tnO  +  tnl  +  next(l  .JttnO)  ♦  tn2  + 

next(2,*tn0)  +  next(2,*tnl)  +  next(2,Jttemp))/8.f ; 
temp  =  next(l ,ttnO) ; 

amuyz2  =  (emu  +  tnO  +  tnl  +  next(l,*tn0)  +  tp2  + 

prev(2, JttnO)  +  prev(2,*tnl)  +  prev(2,Jttemp))/8.f ; 
temp  =  prev(l .JttnO) ; 

amuyz3  =  (emu  +  tnO  +  tpl  +  prev(l,  JttnO)  +  tn2  ♦ 

next(2,*tn0)  +  next(2,Jttpl)  +  next(2,*temp))/8.f ; 
temp  =  prev(l,  JttnO); 

amuyz4  =  (emu  +  tnO  +  tpl  +  prev(l  .JttnO)  +  tp2  + 

prev(2,fttn0)  +  prev(2,Jttpl)  +  prev(2,*temp))/8.f ; 
temp  ■  next  ( 1 ,  JttpO) ; 

amuyzS  =  (emu  +  tpO  +  tnl  +  next(l.ttpO)  +  tn2  + 

next (2,  JttpO)  +  next(2,*tnl)  +  next(2,Jttemp))/8.f ; 
temp  *  next  (1 ,  JttpO)  ; 

amuyz6  =  (emu  +  tpO  +  tnl  +  next (1, JttpO)  ♦  tp2  + 

prev(2, JttpO)  +  prev(2,fctnl)  +  prev(2,Jttemp))/8.f ; 
temp  =  prev(l.ttpO) ; 

amuyz7  =  (emu  +  tpO  +  tpl  +  prev(l  .JttpO)  ♦  tn2  + 

next(2, JttpO)  +  next(2,*tpl)  +  next(2,Jttemp))/8.f ; 
temp  =  prev(l.atpO) ; 

amuyz8  *  (emu  ♦  tpO  +  tpl  +  prev(l ,  JttpO)  +  tp2  + 

prev(2,fctp0)  +  prev(2 ,  Jttpl)  +  prev(2,*temp))/8.f; 

/*  next,  get  the  lambda  value  for  txx,tyy,tzz  positions  */ 
abigul  =  (amuxzl  +  amuxz2  ♦  amuxz3  ♦  amuxz4  + 

amuxyl  +  amuxy2  +  amuxy3  +  amuxy4  + 

amuyzl  +  amuyz2  ♦  amuyz3  +  amuyz4)/12.f ; 
abigu2  =  (amuxzl  +  amuxz3  ♦  amuxzB  +  amuxz6  + 

amuxyl  +  amuxy3  +  amuxyS  +  amuxy6  + 

amuyzS  +  amuyzS  +  amu yz7  ♦  amuyz8)/12.f ; 

alambl  *  alp2ul  -  2.f  *  abigul; 

alamb2  =  alp2u2  -  2.f  *  abigu2; 

txxdx  =  (alp2ul  *  adisp[0]  -  alp2u2  *  adispCl])  + 

(alambl  *  adisp[2]  -  alamb2  *  adisp[3))  + 

(alambl  *  adisp[4]  -  alamb2  *  adisp[5]); 

/*  nos  include  density  term  */ 
txxdx  -  txxdx/rho; 


/*  ==============================================================  */ 

/*  */ 
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/*  Function  calc_txydx 

/* 

void  nodal: :calc_txydx( void) 

{ 

float  anuxyl , amuxy2; 
float  rhov; 
float  tnO.tnl; 
float  adisp[14] ; 

/*  first  gat  tha  displacement  difference  terns  */ 
calc.adisp(adisp) ; 

/*  set  up  temporary  values  */ 
tnO  =  next(0,kcnu) ; 
tnl  =  next(l ,fccmu) ; 

/*  now,  the  nu  for  the  txy  positions  */ 
anuxyl  -  (emu  ♦  next ( 1 ,kcmu))/2 . f ; 
anuxy2  =  (next(O.tcnu)  +  next(l,*tn0))/2.f ; 

txydx  «=  (anuxy2  *  adispClOj  -  anuxyl  *  adisptll])  + 
(anuxy2  *  adisp[12]  -  anuxyl  *  adisp[13]); 

/*  set  up  tenporary  values  */ 
tnO  *  next(O.trho); 
tnl  »  next(l.hrho); 

/*  density  tern  for  y-displacenent  */ 
rhov  ■  (rho  +  tnO  +  tnl  +  next(l,*tn0))/4.f : 

/*  now  include  density  ten  */ 
txydx  =  txydx/rhov; 


/*  SSSSS8SS3S8£SSISSSSSCt«S883SS8C8SI8SC<BCt8tISSS«CSB«SSS38SSC8S  */ 


/*  */ 

/*  Function  calc.txydy  */ 

/*  */ 

/*  . . .  */ 


void  nodal: : cal c_txydy( void) 

float  anuxyl .anuxy 3; 
float  bdisp[14] ; 

/*  first  get  the  displacenent  difference  terns  */ 
calc_bdisp(bdisp) ; 
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/*  now,  th«  bu  lor  the  txy  positions  */ 
aauxyl  =  (can  +  next(l ,Jtcau))/2.f ; 
aauxy3  =  (emu  +  prev(l,fccau))/2.f ; 

txydy  *  (aauxyl  *  bdisp[6]  -  aauxy3  *  bdisp[7])  + 
(aauxyl  *  bdisp[8]  -  aauxy3  *  bdisp[9])  ; 

/*  now  include  density  tera  */ 
txydy  =  txydy/rho; 


/*  */ 

/*  Function  c&lc.txzdx  */ 

/*  */ 


void  model: :calc_txzdx(void) 

{ 

float  aauxzl ,aauxz2; 
float  rhow; 
float  tn0,tn2; 
float  adisp[14] ; 

/*  first  get  the  displaceaent  difference  terms  */ 
calc.adisp(adisp) ; 

/*  set  up  temporary  values  */ 
tnO  *  next(O.hcmu) ; 
tn2  *  next(2,hcmu) ; 

/*  next,  the  mu  value  for  the  txz  positions  */ 

aauxzl  =  (emu  +  tn2)/2.f; 

amuxz2  =  (next(0 , ft emu)  ♦  next(2,*tn0))/2.f ; 

txzdx  *  (amuxz2  *  adisp[6]  -  aauxzl  *  adisp[7]>  + 
(anuxz2  *  adisp[8]  -  aauxzl  *  adisp[9] ) ; 

/*  set  up  temporary  values  */ 
tnO  *  next(O.hrho) ; 
tn2  =  next(2,*rho) ; 

/*  density  tera  for  z-displaceaent  */ 
rhow  =  (rho  ♦  tnO  ♦  tn2  +  next(2,*tn0))/4.f ; 

/*  now  include  density  term  */ 
txzdx  *  txzdx/rhow; 

> 
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*/ 


/*  */ 

/*  Function  calc.txzdz  */ 

/+  */ 

void  nodal : : caic_txzdz(void) 

{ 

float  aauxzl ,aauxz3; 
float  tn2; 
float  edisp [14]; 


/*  first  gat  tha  displacaaant  diffaranca  tarns  */ 
calc.cdisp(cdisp) ; 

/*  sat  up  tamporary  valuas  */ 
tn2  *  naxt(2,ftcau) ; 

/*  naxt,  tha  nu  valua  for  tha  txz  positions  */ 

aauxzl  *  (enu  *  tn2)/2.f; 

aauxz3  1  (emu  +  prav(2,fceau))/2.f ; 

txzdz  *  (aauxzl  *  edisp [6]  -  aauxz3  *  edisp [7] )  + 
(aauxzl  *  cdisp[8]  -  aanxz3  *  edisp [9] ) ; 

/*  now  includa  dansity  ttrs  */ 
txzdz  *  txzdz/rho; 


> 

/*  */ 

/*  Function  calc.tyydy  */ 

/*  */ 

/*  :ic»»tn»nii»untiHi>n«ninimn>uu»<:»i::»:stt  */ 


void  nodal: : calc_tyydy(void) 

float  alp2ul , alp2u3 ; 
float  anuxzl,anuxz2,anuxz3>aauxz4, 
anuxz7 , anuxz8 , aauxz9 , anuxza; 
float  anuxy 1 , anuxy2 , anuxy 3 , anuxy4 , anuxy 7 , anuxy8 ; 
float  anuyzl , anuyz2 , anuyz3 , aauyz4 , anuyzfl , anuyza ; 
float  abigul,abigu3; 
float  alanbl,alanb3; 
float  rhov; 

float  taap,taap2,tsnp3,tsap4; 
float  tnO,tnlltn2,tpO,tpl,tp2; 
float  bdisp[14] ; 

/*  first  gat  tha  displacaaant  diffaranca  tarns  */ 
calc.bdisp(bdisp) ; 
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/*  set  up  temporary  values  */ 
tnO  =  next(0,fcclp2mu) ; 

/*  find  the  lambda+2mu  values  lor  the  tzz,tyy,tzz  positions  */ 
alp2ul  =  (clp2mu  +  tn0)/2.1; 

alp2u3  =  (next (1  ,*clp2mu)  +  next(l , fttnO ))/2.f; 

/*  set  up  temporary  values  */ 

tnO  =  next (0 .ftcmu) ; 

tnl  *  next(l .ftcmu) ; 

tn2  =  next (2, ftcmu) ; 

tpO  =  prev(0 , ft emu) ; 

tpl  =  prev(l.ftcmu) ; 

tp2  =  prev(2 t ft emu) ; 

/*  next,  the  mu  value  lor  the  txz  positions  */ 

amuxzl  =  (emu  +  tn2)/2.1; 

amuxz2  =  (next(0,ftcmu)  ♦  next(2,fttn0))/2.1; 

amuxz3  =  (emu  +  prev(2,ftcmu))/2.1; 

amuxz4  =  (next (0 ,ftcmu)  +  prev(2,fttn0))/2.f ; 

amuxz7  =  (next(l , ft emu)  +  next(2,fttnl))/2.1; 

temp  =  next (1, fttnO) ; 

amuxz8  =  (next(l.fttnO)  +  next(2,fttemp))/2.1; 
amuxz9  =  (next(l,ftcmu)  +  prev(2,fttnl))/2.1; 
temp  =  next(l .fttnO) ; 

amuxza  =  (next(l,fttnO)  +  prev(2,fttemp))/2.1; 

/*  now,  the  mu  for  the  txy  positions  */ 
amuxyl  =  (emu  +  next(l,ftcmu))/2.1; 
amuxy2  ■  (next (0, ftcmu)  +  next (l, fttnO) )/2.f; 
amuxy3  =  (emu  +  prev(l,ftcmu))/2.f ; 
amuxy4  =  (next(0,ftemu)  +  prev(l ,fttn0))/2.1; 
amuxy7  =  (next(i, ftcmu)  ♦  next(l,fttnl))/2.f ; 
temp  =  next (1, fttnO) ; 

amuxy8  *  (next(l ,fttn0)  +  next(l ,*temp))/2.f ; 

/*  finally,  the  mu  for  the  tyz  positions  */ 
temp  =  next(l .fttnO) ; 

amuyzl  =  (emu  +  tnO  +  tnl  +  next (1, fttnO)  +  tn2  ♦ 

next (2, fttnO)  +  next(2,fttnl)  +  next(2,fttemp))/8.1; 
temp  *  next(l, fttnO); 

amuyz2  ■  (emu  +  tnO  +  tnl  +  next(l, fttnO)  +  tp2  + 

prev(2, fttnO)  +  prev(2,fttnl)  +  prev(2,fttemp))/8.1; 
temp  *  pr ev(l, fttnO ); 

amuyz3  *  (emu  +  tnO  +  tpl  +  prev(l, fttnO)  +  tn2  + 

next (2, fttnO)  +  next(2,fttpl)  +  next(2,fttemp) )/8 .1 ; 
temp  *  prev(l, fttnO); 

amuyz4  ■  (emu  +  tnO  +  tpl  +  prev(l, fttnO)  +  tp2  + 

prev(2, fttnO)  ♦  prev(2,*tpl)  +  prev(2,*temp))/8.f ; 
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taap  =  naxt(l ,*tnO) ; 
temp 2  =  naxt(l ,*tn2) ; 
teap3  =  next(l.ktnO); 
tamp3  *  naxt(l  ,ktamp3)  ; 
temp4  =  next(l.ktnO); 

amuyz9  =  (emu  4  next (1 .ktemp)  4  tnl  4  tamp  4  next(l ,ktemp2)  4 
next (2 ,*temp3)  4  next(2,ktnl)  4  next(2,Jttemp4))/8.1; 
tamp  =  next(l.ktnO); 
temp 2  -  next(l.ktnl); 
tamp3  =  naxtd.ktnO); 
tamp3  =  next (1 ,ttemp3) ; 
temp4  =  naxtd.ktnO); 

amuyza  =  (naxt(l.ktnl)  4  next(l.ktemp)  4  tnl  4  naxt(l.ktnO)  4 
prev(2 ,fctemp2)  4  prav(2,fctemp3)  4  prev(2,ktnl)  4 
prev(2,*temp4))/8.1 ; 

/*  naxt,  gat  ttaa  lambda  value  lor  txx.tyy.tzz  positions  */ 
abigul  =  (amuxzl  4  amuxz2  4  amuxz3  4  amuxz4  4 

amuxyl  4  amuxy2  4  amuxy3  4  amuxy4  4 

amuyzl  4  amuyz2  4  amuyz3  4  amuyz4)/12.f ; 
abigu3  -  (amuxz7  4  amuxz8  4  amuxzB  4  amuxza  4 

amuxy7  4  amuxy8  4  amuxyl  4  amuxy2  4 

amuyzl  4  amuyz2  4  amuyzS  4  amuyza)/12.1; 

alambl  =  alp2ul  -  2.1  *  abigul; 

alamb3  *  alp2u3  -  2.1  *  abigu3; 

tyydy  *  (alamb3  *  bdisptO]  -  alambl  *  bdispCl])  4 

(alp2u3  *  bdisp[2]  -  alp2ul  a  bdi«p[3])  4 

(alamb3  *  bdisp[4]  -  alambl  *  bdisp[B] ) ; 

/*  aat  up  temporary  value*  */ 
tnO  =  next (O.fcrho) ; 
tnl  *  next(l.trho); 

/*  density  term  lor  y-displacement  */ 
rhov  *  (rbo  4  tnO  4  tnl  4  next(l,ktn0))/4.1 ; 

/*  now  include  density  term  */ 
tyydy  *  tyydy/rhov; 


/*  SSISSCSXSSSSSSISS8CS88tntS£RB|C88lttltll8Stl8»28KSSSISSSCSC  */ 


/*  */ 

/*  Function  calc.tyzdy  */ 

/*  */ 

/*  u***mm**mmm*mmm********nmmm**u*mnmmmmmmmmm**n*mBmm*xn*mm*m*m  */ 

void  model: :calc_tyzdy(void) 
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< 

float  amuyzl ,amuyz3; 
float  rhow; 
float  temp; 

float  tnO.tnl ,tn2,tpl ; 
float  bdisp[14] ; 

/*  first  gat  the  displacement  difference  terms  */ 
calc_bdisp(bdisp) ; 

/*  set  up  temporary  values  */ 

tnO  -  next(0,kcmu) ; 

tnl  =  next(l ,fccmu) ; 

tn2  =  next (2 ,fccmu) ; 

tpl  =  prev(l ,kcmu) ; 

/*  finally,  the  mu  for  the  tyz  positions  */ 
temp  -  next ( 1 , fctnO) ; 

amuyzl  =  (emu  +  tnO  +  tnl  +  next(l.ktnO)  +  tn2  + 

next(2,*tn0)  +  next(2,*tnl)  +  next(2,ktemp))/8.f ; 
temp  -  prev(l.ktnO) ; 

amuyz3  =  (emu  +  tnO  +  tpl  +  prev(l.ktnO)  +  tn2  + 
next (2 ,ktnO)  +  next(2,ktpl)  +  next(2,ktemp))/8.f ; 

tyzdy  x  (amuyzl  *  bdispClO]  -  amuyz3  *  bdisptll] )  ♦ 

(amuyzl  *  bdisp[12]  -  amuyz3  *  bdisp[l3]); 

/*  set  up  temporary  values  */ 
tnO  x  next(O.krho) ; 
tn2  =  next  (2,  Jtrho) ; 

/*  density  term  for  z-displacement  */ 
rhow  «  (rho  +  tnO  +  tn2  +  next (2,ktn0))/4.f ; 


/*  now  include  density  term  */ 
tyzdy  =  tyzdy/rhow; 

> 


/* 

*/ 

/* 

Function  calc.tyzdz 

*/ 

/* 

*/ 

void  model: :calc_tyzdz(void) 
■C 

float  amuyzl ,amuyz2; 
float  rhov; 
float  temp; 

float  tn0,tnl,tn2,tp2; 
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float  cdisp[14]; 

/*  first  get  the  displacement  difference  terms  */ 
calc.cdisp(cdisp) ; 

/*  set  ap  temporary  values  */ 

tnO  =  next (0, ft emu) ; 

tnl  =  next ( 1 , ftcmu) ; 

tn2  =  next (2, ft emu) ; 

tp2  =  prev(2,ftcmu) ; 

/*  finally,  the  mu  for  the  tyz  positions  */ 
temp  *  next(l ,fttnO) ; 

amuyzl  =  (emu  +  tnO  +  tnl  +  nextd, fttnO)  +  tn2  + 

next(2,fttn0)  *  next(2,fttnl)  +  next (2,fttemp))/8.f ; 
temp  -  next(l .fttnO) ; 

amuyz2  =  (emu  +  tnO  +  tnl  +  next (1, fttnO)  +  tp2  + 

prev(2, fttnO)  +  prev(2,fttnl)  +  prev(2,fttemp))/8.f ; 

tyzdz  =  (amuyzl  *  cdispClO]  -  amuyz2  *  cdispUll])  + 

(amuyzl  *  cdisp[l2]  -  amuyz2  *  cdisp[13]); 

/*  set  up  temporary  values  */ 
tnO  =  next(O.ftrho) ; 
tnl  *  next(l,ftrho); 

/*  density  term  for  y-displacement  */ 
rhov  *  (rho  +  tnO  +  tnl  +  next (1, fttnO) )/4.f; 

/*  now  include  density  term  */ 


tyzdz  *  tyzdz/rhov; 

> 

/ *  csnmntnsusnsRinciMicncninncKcicnnnKsncsn  * / 

/*  */ 

/*  Function  calc.tzzdz  */ 

/*  */ 


void  model: : calc.tzzdz (void) 

{ 

float  alp2ul,alp2u4; 

float  amuxzl ,amuxz2>amuxz3,amuxz4,amuxzb,amuxzc; 
float  amuxyl,amuxy2,aatuxy3lamuxy4> 
amuxyS , amuxya , amuxyb , amuxy c ; 
float  amuyz 1 , amuyz2 , amuyz3 , amuyz4 , amuyzb , amuyze ; 
float  abigul,abigu4; 
float  alambl,alamb4; 
float  rhow; 

float  temp,temp2,temp3,temp4; 
float  tn0,tnlltn2,tpl,tp2; 
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float  cdisp[14] ; 


/*  first  get  the  displacement  difference  terms  */ 
calc.cdisp(cdisp) ; 

/♦  set  up  temporary  values  */ 
tnO  =  next(0,kclp2mu) ; 

/*  find  the  lambda+ 2mu  values  for  the  txx.tyy.tzz  positions  */ 

alp2ul  =  (clp2mu  +  tn0)/2.f; 

alp2u4  =  (next(2,kclp2mu)  +  next(2,ktn0))/2.f ; 

/*  set  up  temporary  values  */ 

tnO  =  next (0 ,kcmu) ; 

tnl  =  next(l ,kcmu) ; 

tn2  -  next(2,kcmu) ; 

tpl  =  prev(l ,kcmu) ; 

tp2  =  prev(2,kcmu) ; 

/*  next,  the  mu  value  for  the  txz  positions  */ 
amuxzl  »  (emu  +  tn2)/2.f; 

amuxz2  =  (next(O.kcmu)  +  next(2,ktn0))/2.f ; 
amuxzS  *  (emu  +  prev(2,kcmu))/2.f ; 
amuxz4  *  (next(O.kcmu)  +  prev(2,ktn0))/2.f ; 
amuxzb  =  (next (2,kcmu)  +  next(2,ktn2))/2.f ; 
temp 2  *  next (2 ,ktnO) ; 

amuxzc  *  (next(2,ktn0)  +  next(2,ktemp))/2.f ; 

/*  non,  the  mu  for  the  txy  positions  */ 
amuxyl  *  (emu  +  next(l ,kcmu))/2.f ; 
amuxy2  =  (next(0,fccmu)  +  next(l,ktn0))/2.f ; 
amuxy3  =  (emu  +  prev(l ,kcmu))/2.f ; 
amuxy4  =  (next(0,kemu)  +  prev(l ,fctn0))/2.f ; 
amuxyS  *  (next(2,kcmu)  +  next(2,ktnl))/2.f ; 
temp  -  next(i,ktnO); 

amuxya  »  (next(2,ktn0)  +  next(2,ktemp))/2.f; 
amuxyb  =  (next(2,kcmu)  +  next(2,ktpl))/2.f ; 
temp  ■  prev(l,ktnO); 

amuxyc  =  (next(2,ktn0)  +  next(2,ktemp))/2.f ; 

/*  finally,  the  mu  for  the  tyz  positions  */ 
temp  *  next(l,ktnO); 

amuyzl  »  (emu  ♦  tnO  ♦  tnl  +  nextd.ktnO)  +  tn2  + 

next(2,ktn0)  ♦  next(2,ktnl)  +  aext(2,ktemp))/8.f ; 
temp  *  next(l.ktnO); 

amuyz2  *  (emu  +  tnO  +  tnl  +  next(l,ktnO)  +  tp2  + 

prev(2,ktn0)  +  prev(2,ktnl)  +  prev(2,ktemp))/8.f ; 
temp  *  prev(l,ktnO); 

amuyz3  *  (emu  ♦  tnO  +  tpl  +  prev(l.ktnO)  +  tn2  + 
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next(2,fctn0)  +  next(2,fctpl)  +  next (2 ,*temp) )/8 .1 ; 
temp  =  prev(l ,fctnO) ; 

amuyz4  =  (can  +  tnO  +  tpl  +  prev(l.htnO)  +  tp2  + 

prev(2,*tn0)  +  prev(2,*tpl)  +  prev(2,*teap))/8.f ; 
temp  =  next(2,*tnO) ; 
teap2  =  next(2,*tnl) ; 
teap3  =  next(l ,*tnO) ; 
temp3  =  next(2,fctemp3) ; 
teap4  =  next(l ,*tnO) ; 

amuyzb  =  (next(2,*tn2)  +  next(2,Jrtemp)  +  next (2 ,ktemp2)  + 

next (2,ktemp3)  +  tn2  +  next(2,*tn0)  +  next(2,fctnl)  + 
next(2,fctemp4))/8 .f ; 
temp  s  next(2,*tn0) ; 
temp2  3  next(2,fctpl) ; 
temp3  3  pr*v(l,*tnO); 
temp 3  >  next(2,kteap3) ; 
temp4  3  prevd.fttnO); 

amuyzc  *  (next (2 ,*tn2)  ♦  next(2 ,*temp)  +  next(2,*temp2)  + 

next(2,fctemp3)  +  tn2  ♦  next(2,*tn0)  +  next(2,fctpl)  + 
next(2,ttemp4))/8.f ; 

/*  next,  get  the  lambda  value  for  txx.tyy ,tzz  position*  */ 
abigul  *  (amuxzl  ♦  amnxz2  +  aauxz3  +  aantxz4  ♦ 

amnxyl  +  aauxy2  +  aanxy3  +  aauxy4  + 

aauyzl  +  amnyz2  +  amnyz3  +  asruyz4)/12.f ; 
abigu4  3  (amuxzl  +  amuxz2  ♦  amuxzb  +  amuxzc  + 

amuxy9  +  amuxya  +  aauxyb  +  amuxyc  + 

aauyzl  +  aauyz3  +  amuyzb  +  amuyzc) /12.f ; 

alaabl  3  alp2ul  -  2.1  *  abigul; 

alamb4  3  alp2u4  -  2.1  *  abigu4; 

tzzdz  3  (alaab4  *  cdispCO]  -  alambl  *  cdisptl])  + 

(alamb4  *  cdisp[2]  -  alambl  *  cdisp[3])  ♦ 

(alp2u4  *  cdisp[4]  -  alp2ul  *  cdisp[5])  ; 

/*  set  up  temporary  values  */ 
tnO  3  next(O.trho) ; 
tn2  3  next(2,*rho) ; 

/*  density  term  for  z-displacement  */ 
rhow  *  (rho  +  tnO  +  tn2  +  next(2,ktn0))/4.f ; 

/*  now  include  density  term  */ 
tzzdz  3  tzzdz/rhow; 


/♦ 


*/ 
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/*  */ 

/*  Function  divcurl  -  calculate  the  divergence  and  curl  */ 

/*  +/ 

void  aodel: : divcurl (void) 

{ 

int  ix.iy.iz; 

mono  float  div_ain,div_aax; 

/*  zero  out  div  and  curl  values  */ 
div  =  curlx  *  curly  =  curlz  =  0.0; 

/*  calculate  divergence  and  curl  */ 
ix  =  this.coordinate(O) ; 
iy  *  this_coordinate(l) ; 
iz  =  this_coordinate(2) ; 

if ((ix  <  1  II  ix  >=  IDIM)  II  (iy  <  1  II  iy  >=  YDIM)  II 
(iz  <  1  II  iz  >-  ZDIM) ) 


else 

div  =  (next(O.fcul)  -  ul)/dx  +  (vl  -  prev(l ,*vl))/dy  + 

(wl  -  prev(2,kvl))/dz; 

curlx  =  (next ( 1 ,kwl)  -  wl)/dy  -  (next(2,Jtvl)  -  vl)/dz; 
curly  =  (next (2, Jrul)  -  ul)/dz  -  (»i  -  prev(0 ,kwl) )/dx; 
curlz  *  (vl  -  prev(0,kvl))/dx  -  (next(l.kul)  -  ul)/dy; 

> 

/*  get  ainiaua  and  naxiatun  div  */ 
div.min  <?=  div; 
div.aax  >?=  div; 

printf ("ainiaua  div:  %14.7e  aaxiaua  div:  %14.7e\n",  div_ain,div_aax) ; 

> 

/*  =s::s:s»:ss:ss:sxs3s:ns:c3»3»csin«»nn»tM»:sKsss:3  */ 

/*  */ 

/*  Function  force  */ 

/*  */ 

/*  ss::::3:cnsscz«scsU3>it«KnsnentHn»>n>s«3H»3>in«:  */ 

force(int  ix , int  jy,int  kz.int  ltiae, float  efsorsu, float  efsorsv, float  efsorsw) 

C 

int  del; 

float  rn,rv,rv,hx,hy,hz; 

float  xpvidt , gtp , gt , dircsx , dircsy , dir csz ; 

float  tiae  *  ltiae  *  dt; 

float  hh( float  r); 

del  *  6  *  dx; 
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/ *  calculate  range  values  (radii)  for  u,v,  and  w  grid  points  */ 
ru  =  sqrt((ix*ix)  *  (dx*dx)  +  (jy*jy)  *  (dy*dy)  +  (kz*kz)  *  (dz*dz)); 
rv  =  sqrt(((ix+0.5f )*(ix+0.6f ))  *  (dx*dx)  +  ((jy+0.6f)*(jy+0.6f))  * 
(dy*dy)  +  (kz*kz)  *  (dz*dz)); 

rv  *  sqrt ( ( ( ix+O . 6f )*(ix+0.5f ))  *  (dx*dx)  +  (jy*jy)  *  (dy*dy)  + 

( (kz+O .Si)* Ocz+0 . 5f ) )  *  (dz*dz)) ; 

/*  calculate  ’jagged*  sine  save  for  x,y,  and  z  -  spatial  scaling  function 
applied  to  the  tine  dependent  source  function  */ 
hx  =  2.f/(del*del)  *  (-ru*hh(ru)  +  2.f  *  (ru  -  del/2. f )*hh(ru-del/2.f )  - 
(ru-del)ehh(ru-del)) ; 

hy  «  2.f/(del*del)  •  (-rv*hh(rv)  +  2.f  *  (rv  -  del/2. f )*hh(rv-del/2.f )  - 
(rv-del)*hh(rv-del)) ; 

hz  *  2.f/(del*del)  *  (-rv*hh(rv)  +  2.f  *  (rv  -  del/2. f)*hh(rv-del/2.f)  - 
(rv-del)*hh(rv-del)) ; 

/*  now  the  tine  dependence  (gaussian)  */ 
xpvidt  =  (PI  *  fpeak)  *  (PI  *  fpeak); 
gtp  =  exp(  -  (time* tine)  *  xpvidt); 
gt  *  -2 .f *xpvidt*(l .f-2 .f *xpvidt*(tine*tine) )  *  gtp; 

/*  finally,  put  it  all  together  */ 
dir csx  =  (ix*dx)/ru; 
dircsy  =  ((jy+0.5f )*dy)/rv; 
dircsz  =  ((kz+0.5f)*dz)/rv; 

sfsorsu  s  gt  *  hx  *  dir csx; 

♦fsorsv  *  gt  *  hy  *  dircsy; 
efsorsv  =  gt  *  hz  *  dircsz; 

> 

/*  SSS::SSSSSSSSSCSSS88S&SCSSSS8«CSSXXSSS8C»acSBSSSSeSSS 

/* 

/*  Function  hetiso.cs 

/* 

/*  :«s:i»tn>isunue»i:»sscein»iisiii»sst:E::ts> 

void  hetisoO 

nono  int  It ins; 
nono  char  lsors; 

void  source(int  ltine); 
void  output .files (nono  int  ltine); 

lsors  ■  TRUE;  /*  input  a  source  •/ 

/*  boundary  conditions  can  be  set  up  outside  tine  loop  for  parallel  */ 

/*  processing  -  set  up  for  asorbing  boundaries  */ 
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[domain  modal] . {boundary  C ) ; > 


/*  time  loop  -  start  at  time  step  nti  =  3  */ 
lor  (ltime  =  ntl;  ltime  <=  nt;  ltime++) 

{ 

/*  calculate  the  stress  derivative  terms  */ 
[domain  model] .{stress () ;> 

/*  now,  calculate  displacements  */ 

[domain  model] . {calc_displ ( ) ; > 

/*  update  the  time  series  */ 

[domain  model] .{update();> 

/*  input  the  source  */ 
il(lsors) 

source (ltime) ; 

/*  output  options  here  */ 
il  (ltime  '/.  snap.file  ==  0) 
output .files (ltime) ; 

/*  output  to  frame  buffer  */ 
if  (ltime  */.  snap.fb  «  0) 

[domain  model] .{output_fb() ;} 

>  /*  end  of  time  loop  */ 

> 


/*  ::»:::cs<s>«nnuni>cau»»n»»n»»H:»st»UK»st  */ 
/*  */ 
/*  Function  hh  */ 
/*  */ 
/*  SIBSSSBBICSBSUESttsnUSSBtnCSSlBSBtttBBSIStBSiacISSSSBtBBSS  $  / 


float  hh(float  r) 

{ 

if (r  <  0.0) 

return(0.0f ) ; 

else 

{ 

if (r  <  0.00001) 
return (0.6) ; 

else 

retumd  .0) ; 

> 

> 

/*  •/ 
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*/ 

*/ 

♦/ 


/*  Function  initialize  -  initialize  the  grid 

/* 

void  model: : initializeO 

{ 

uO  =  ul  =  u2  =  0.0; 
vO  =  vl  =  v2  =  0.0; 
wO  =  ml  ~  h2  =  0.0; 

> 


/*  */ 

/*  Function  input.defs  *■/ 

/*  +/ 

void  input _defs() 

FILE  *fp; 

int  nz,i, j ,k,pr_indx; 


int  nzonee , nil ,nyl ,nzl ,nxn,nyn,nzn; 
float  dt2,cclamb,ccmu,trho; 

/*  open  the  file  containing  model  parameters  */ 
if  ((fp  =  fopenCxmodel.dat", "r")}  **  *ULL) 

printf ("Error  opening  input  file:  model. dat\n"); 

/*  read  model  parameters  */ 
fscanf (fp,"%d",  tat); 

printf (“lumber  of  time  steps:  %d\n\n" ,  nt); 
f scanf (fp,"XfWfXf" ,  *dx ,*dy ,*dz,*dt) ; 
fscanf (fp,"%d%d",  tanap.file,  tanap.fb) ; 
f scanf (fp, a'XdXd%d%d" ,  taut _fb, taut _xy , tout _yz, tout _xz) ; 
f scanf  (f p , "W/Md'At" ,  hisx.hisy ,*isz,*fpeak) ; 

ntl  *  3; 

/*  calculate  lambda  factors  for  use  in  ddispl  */ 

dt2  *  dt*dt; 

xlamb  -  dt2/(dx*dx) ; 

ylamb  «  dt2/(dy*dy); 

zlamb  =  dt2/(dz*dz); 

xylamb  *  dt2/(dx*dy); 

xzlamb  *  dt2/(dx*dz) ; 

yzlaab  *  dt2/(dy*dz) ; 

/*  read  in  number  of  zones  of  different  elasticity,  then  read  in 
constants  for  each  zone  */ 
fscanf (fp,”%d",  tazones); 

for(nz*0;  nz<nzones;  nz++) 
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{ 

f scant (fp, "V.d'Al'/.d’/.dy.d'/.d" ,  Anzl ,tayl ,*nzl , knxn,*nyn,*nzn) ; 
f scant (fp,  "’/.f'/.f'/.f  " ,  Acclamb,ftccmu,Atrho) ; 

/*  make  the  parallel  variable  assignments  in  parallel  */ 
[domain  model] . { 
int  iz.iy.iz; 

iz  =  this.coordinate(O) ; 
iy  =  this.coordinate(l) ; 
iz  =  this_coordinate(2) ; 

/*  select  the  active  processors  */ 
il((iz  >=  ail)  kk  (iz  <=  nzn)  U 
(iy  >=  nyi)  A*  (iy  <=  nyn)  Aft 
(iz  >=  nzl)  kk  (iz  <=  nzn)) 

clamb  »  cclamb ; 
emu  *  ccmu ; 

dp2mu  =  cclamb  +  2.f  *  ccmu; 
rho  *  trho; 

>  /*  end  it  */ 

>  /*  end  domain  */ 

>  /*  end  for  */ 

> 


/*  ■■■■■■■■■•■•■•■■■■■■■■■■■■■••••■■■••■•■■■■■■■■•••■■■■■■•■■■■■I*  */ 
/*  */ 

/*  Function  itoa  -  integer  to  ASCII  conversion  */ 

/*  */ 

/*  »»»tnainnnnnmmanuinaninBnmmiamnai  */ 

void  itoa(int  n,  char  sD) 

i 

int  i,  sign; 

void  reverse(char  s  □  ) ; 


if ((sign  *  n)  <  0)  /*  record  sign  */ 

n  =  -n;  /*  make  n  positive  */ 

i  *  0; 
do  { 

s[i++]  *  n  %  10  +  *0*; 

>  uhile  ((n  /-  10)  >  0); 
if  (sign  <  0) 

• [i++]  * 
s  Ci]  *  *\0' ; 
reverse(s) ; 

> 

/*  «<nunHi>nm»imiMiMitiini»nMniHini>n«nn>tt  */ 

54 


/*  */ 

/*  Function  output_fb  -  output  to  fraae  buffer  */ 

/*  */ 

/*  output  options:  */ 

/*  out.fb  =  0  for  no  fraae  buffer  output  */ 

/*  *  1  to  display  xy  plane  at  ZDIH/2-1  */ 

/*  =  2  to  display  yz  plane  at  XDIH/2-1  */ 

/*  s  3  to  display  xz  plane  at  YDIM/2-1  */ 

/*  */ 

/  *  :SSSSSSSSKS8SSBC»8B8S8se:SSSSSSSSSSSB8SIISS:SCSCX8Si:SS»=:5  */ 

void  model: : output _fb(void) 


mono  int  nx2  =  (XDIM/2)-l,  ny2  =  (YDIK/2)-l,  nz2  =  (ZDlM/2)-l; 
mono  float  slope,  intercept; 
unsigned  char  color; 
int  ix,iy,iz; 

/*  calculate  divergence  and  curl  */ 
divcurlO ; 

/*  get  ainiaua  and  aaxiaua  div  and  curl  via  reduction,  use  to  */ 
/*  scale  for  fraae  buffer  */ 

slope  *  23. Of  /  ((>?«  div)  -  (<?*  div)); 
intercept  ■  l.Of  -  elope  *  (<?*  div); 

/*  calculate  the  color  */ 
color  *  slope  *  div  *  intercept; 

/*  select  processors  to  display  */ 
ix  *  this.coordinate(O); 
iy  =  this.coordinate(l) ; 
iz  =  this_coordinate(2) ; 

/*  plot  the  xy  plane  */ 
if(out_fb  **  1) 

< 

if(iz  **  nz2) 

plot_x_y( (unsigned  short)ix,  (unsigned  short)iy,  color); 

> 

/*  plot  the  yz  plane  */ 
if(out_fb  ■*  2) 

{ 

if(ix  **  nx2) 

plot_x_y( (unsigned  short)iy,  (unsigned  short)iz,  color); 

> 

/*  plot  the  xz  plane  */ 
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if (out_fb  ==  3) 

{ 

if (iy  ==  ny2) 

plot_x_y( (unsigned  short )ix,  (unsigned  short)iz,  color); 

> 


/*  */ 

/*  Function  output _files  -  write  out  ASCII  diak  files  */ 

/*  */ 

/*  output  options:  ♦/ 

/*  xy_out  *  1  to  output  div_xy  and  curlz  ♦/ 

/♦  yz_out  *  1  to  output  div_yz  and  curlx  ♦/ 

/♦  xz_out  *  1  to  output  div_xz  and  curly  ♦/ 

/*  ♦/ 

/♦  sniuu«s:tttsnc>Ks»narnict:ii>ii>««nHMKcntK»:i::  ♦/ 

void  output_f iles(nono  int  It ins) 

< 

■ono  int  i, j ,k,indx,indx2; 


■ono  int  nx2  *  (XDIM/2)-l,  ny2  =  (YDIM/2)-l,  nz2  =  (ZDIH/2)-l; 
char  divxy [15],  divyz[15],  divxz[15]; 
char  cur lx [15],  curly [16],  curlz[16] ; 
char  tatep[4] ,  extent [5]; 

FILE  *fp_divxz ,  efp_divyz,  *fp_divxy ; 

FILE  efp_curlx,  *fp_curly,  *fp_curlz; 
char  ♦street (char  efilel,  char  efile2); 
char  *etrcpy(char  efilel,  char  efile2); 
void  itoa(int  n,  char  es) ; 

strepy ( extent , " . dat " ) ; 
itoa(ltiaie,  tstep); 

itoa(lti*e,  tstep); 

/*  open  ASCII  output  files  for  div  and  curl  */ 
if (out_xy) 

strepy (divxy,"divxy_") ; 

strcat(divxy,  tstep); 

strcat(divxy,  extent); 

printf ("Opening  file  %s\n”,  divxy); 

if  ((fp.divxy  «  fopen(divxy,"w"))  ■■  MULL) 

printf ("Error  opening  output  file:  %s\n",  divxy); 
strepy (cur lz, "cur lz_") ; 
street (curlz,  tstep); 
strcat(curlz,  extent); 
printf ("Opening  file  %s\n",  curlz); 
if  ((fp.curl z  *  fopen(cuxlz,"w"))  *■  TOLL) 
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print! ("Error  opening  output  file:  %e\n" ,  curlz) ; 

> 

if(out.yz) 

strcpy(divyz, "divyz.") ; 

•treat (divyz ,  tetep) ; 

•treat (divyz,  extent); 

print! ("Opening  file  Xs\n",  divyz); 

if  ((fp.divyz  *  fopen(divyz,"s"))  *«  TOLL) 

print! ("Error  opening  output  file:  X«\n" ,  divyz); 
•trcpy(curlx,"curlx_") ; 
strcat(curlx,  tetep); 

•treat ( cur lx ,  extent ) ; 

print! ( "Opening  file  Xs\n",  curlx); 

if  ((fp_curlx  *  fopen(curlx,"u"))  *■  IULL) 

print! ("Error  opening  output  file:  */.»\n",  curlx); 

> 

if (out_xz) 

{ 

■trcat(divxz, "divxz_") ; 

■treat (divxz,  tetep); 

•treat (divxz,  extent); 

print! ("Opening  file  Xs\n",  divxz); 

if  ((fp_divxz  «  fopen(divxz,“v"))  **  TOLL) 

print! ("Error  opening  output  file:  Xs\n",  divxz); 

■treat (curly , "curly.") ; 

■treat (curly ,  tstep); 

■treat (curly,  extent); 

print! ("Opening  file  Xs\n",  curly); 

if  ((fp.curly  *  f open ( cur ly,"w"))  ■■  TOLL) 

print! ("Error  opening  output  file:  Xs\n",  curly); 

> 

/*  calculate  divergence  and  curl  */ 

[domain  model] . {divcurl ( ) ; } 

/*  write  out  principal  axis  output  plot  files  (snapshot  output  plots)  */ 
if (out.xy) 

fprintf (fp.divxy , "divergence  -  xy  plane  at  z  *  Xd:\n",  nz2); 
for(i  ■  0;  i  <  XDM;  i*+) 
for(j  ■  0;  j  <  TOIM;  j++) 

indx  ■  index_froa.gr id (i,nz2, j) ; 

fprintf (fp.divxy,  "X9.2e  ", point [indx] .div); 

if(((j+l)  X  8)  ■■  0)  fprintf (fp.divxy,  "\n"); 

} 

fprintf (fp.divxy,  "\n"); 
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iprintl (ip.curlz,  "curlz  -  xy  plana  at  z  *  54d:\n",  nz2); 
lor(i  =  0;  i  <  XDIM;  i++) 

lor( j  =  0;  j  <  YDIM;  j++) 

{ 

indx  ®  index_irom_grid(i,nz2,  j)  ; 

iprintl (Ip.curlz,  "*/,9 . 2e  “.point [indx]  .curlz); 

il(((j+l)  54  8)  ==  0)  Iprintldp.curlz,  "\n"); 

> 

iprintl (Ip.curlz ,  "\n"); 

> 

ii (out_yz) 

{ 

Iprintl (Ip.divyz,  "divergence  -  yz  plane  at  x  *  54d:\n",  nx2) ; 
ior(i  =  0;  i  <  YDIM;  i++) 

lor( j  *  0;  j  <  ZDIM;  j++) 

{ 

indx  -  index_irom_grid(i,nx2, j) ; 

Iprintl  (ip_divyz ,  '"/.9 . 2e  ",  point  [indx]  .  div)  ; 
il(((j+l)  54  8)  *=  0)  Iprintl (Ip.divyz,  "\n"); 

> 

Iprintl (Ip.divyz,  "\n"); 

Iprintl (ip_curlx,  "curlx  -  yz  plane  at  x  *  54d:\n",  nx2); 
ior(i  =  0;  i  <  YDIM;  i++) 

ior(j  *  0;  j  <  ZDIM;  j++) 

{ 

indx  =  index_lron_grid(i,nx2, j) ; 

Iprintl (ip_cur lx,  "549. 2e  " .point [indx] . curlx) ; 
i*(((j+l)  54  8)  «  0)  Iprintldp.curlz,  "\n"); 

> 

Iprintl (Ip.curlx ,  "\a”); 

> 

ii(out_xz) 

{ 

Iprintl (ip.divxz,  "divergence  -  xz  plane  at  j  =  54d: \n" ,  ny2) ; 
lor(i  *  0;  i  <  XDIM;  i++) 

lorCj  «=  0;  j  <  ZDIM;  j++) 

{ 

indx  *  index_iro»_grid(i,ny2, j) ; 

Iprintl (Ip.divxz ,  "549. 2e  " .point [indx] .div) ; 
ii(((j+l)  54  8)  ■«  0)  Iprintl (ip.divxz,  "\n"); 

> 

Iprintl (Ip.divxz ,  "\n") ; 

iprintl (ip.curly,  "curly  -  xz  plane  at  y  ■  54d:\n",  ny2); 
lor(i  *  0;  i  <  XDIM;  i++) 


58 


for(j  =  0;  j  <  ZDIM;  j++) 

{ 

indx  =  index_from_grid(i,ny2, j) ; 
fprintf(fp_ curly ,  "'/.9.2e  *' .point [indx]  . curly) ; 
ii(((j+l)  */.  8)  ==  0)  fprintf(fp_ curly,  "\n"); 

> 

fprintf (fp.curly,  "\n”); 

> 

> 


/*  */ 

/*  Function  reverse  */ 

/*  */ 

/*  :=8s:scsscE=zs:s=:ss===sss==sss:=:scsssss»:xssEzcs£s=ssss===s  * / 

void  reverse(char  b □ ) 

{ 

int  c.i.j; 

int  strlenCchar  *s) ; 


lor  (i  *  0,  j  *  strlen(s)-l;  i<j;  i++,  j — ) 

{ 

c  ■  ■  [i]  ; 

•  til  ■  ; 

»[j]  =  c; 

> 

> 

/••iMHnsunnaifiniminniimininiiniMinnniBmn  */ 


/*  */ 

/*  Function  setup.fb  -  initialize  frane  buffer  */ 

/*  and  define  color  table  */ 

/*  */ 


/$SSSttX8ilBSSII«»SI8SSISS8B8SII«IISIStSllt«BtSISSSIB»ISSt«CSSI  */ 

void  setup.fb(void) 

{ 

int  i; 

/*  initialize  the  fraae  buffer  e/ 
init_fra*e_buff er(128, 128) ; 

/*  define  the  color  table  «/ 
for(i»0;  i<26;  i**) 

set_color(i+l,  240-i*10,  i*10,i*10); 

> 

/*  n::i3:a:E»rKE::«snHisH»uu*iHiii»<iHH>3H3t>n«  */ 

/*  */ 

/*  Function  source  -  updates  the  displacement  values  */ 

/*  for  the  source  by  calling  function  force  */ 
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/*  which  calculates  the  force  representation  of  */ 

/*  a  point  source.  The  source  location  is  given  */ 

/*  by  [isx]  [isy]  [isz] .  */ 

/*  */ 

void  source(int  ltime) 

< 

int  i,j,k,ix,jy,kz, pt.indx; 
float  xoff  =  1.0e-10; 
float  f sorsu.f sorsv,! sorsv; 

for  (i  =  0;  i  <=  DEL;  i++) 

ix  =  0  -  i; 

for  (j  *  0;  j  <*  DEL;  j++) 
jy  -  0  -  j; 

for  (k  =  0;  k  <=  DEL;  k++) 
kz  =  0  -  k; 

if  ((ix  «■  0)  kk  (jy  «■  0)  kk  (kz  »*  0)) 

• 

•1  99 

f  or ce ( ix , j  y , kz , ltime ,  kt  sorsu ,  kt  sorsv , ki sorsv) ; 

/*  add  the  forcing  term  to  displacements  at  corners  of  cubes  */ 
if  ((ix  !*  0)  kk  (jy  !«  0)  kk  (kz  !•  0)) 

/*  get  the  point  index  for  the  location  to  be  calculated  */ 
pt.indx  ■  index_fron_grid(isx+ix,isy+jy, isz+kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  +  fsorsu; 
pt.indx  *  index_fron_grid(isx+ix, isy-jy, isz+kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  +  fsorsu; 
pt.indx  »  index_fron_grid(isx+ix, isy-jy, isz-kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  +  fsorsu; 
pt.indx  «  index.from_grid(isx+ix,isy+jy, isz-kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  +  fsorsu; 
pt.indx  «  index.! rom.gr id( isx-ix, isy- jy , isz+kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  -  fsorsu; 
pt.indx  *  index_fron_grid(isx-ix,isy-jy,isz-kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  -  fsorsu; 
pt.indx  *  index_from_grid(isx-ix,isy+jy, isz+kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  -  fsorsu; 
pt.indx  =  index.f rom.gr id( isx-ix , isy+ jy, isz-kz) ; 
point  [pt.indx]  .ul  ■  point  [pt.indx]  .ul  -  fsorsu; 

pt.indx  *  index.f ron.gr id(isx+ix,isy+jy, isz+kz); 
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point [pt.indx] .vl  =  point [pt.indx] . vl  +  fsorsv; 
pt_indx  =  index_fron_grid(isx+ix,isy-jy-l,i8z+kz) ; 
point Cpt_indx] .vl  =  point [pt.indx] .vl  -  fsorsv; 
pt_indx  =  index_fron_grid(isx-ix-l,isy+jy,isz+kz); 
point [pt.indx] .vl  =  point [pt.indx] . vl  +  fsorsv; 
pt.indx  =  index_fron_grid(isx-ix-l,isy-jy-l,isz+kz) ; 
point [pt_indx3 .vl  =  point [pt_indx3 .vl  -  fsorsv; 
pt_indx  =  indsx_froB_grid(isx+ix,isy+jy,isz-kz) ; 
point [pt_indx3 .vl  =  point [pt.indx] .vl  +  fsorsv; 
pt.indx  =  indsx_from_grid(isx+ix,isy-jy-l ,isz-kz) ; 
point [pt_indx3 .vl  =  point [pt_indx3 .vl  -  fsorsv; 
pt.indx  =  ind«x_froB_grid(isx-ix-l,isy+jy ,isz-kz) ; 
point [pt_indx3 .vl  =  point Cpt_indx3 .vl  +  fsorsv; 
pt.indx  =  indsx_fron_grid(isx-ix-i,isy-jy-l,isz-kz) ; 
point [pt_indx3 .vl  =  point [pt_indx3 .vl  -  fsorsv; 

pt.indx  =  index_fron_grid(isx+ix,isy+jy,isz+kz) ; 
pointCpt_indx3.nl  =  pointCpt_indx3.nl  +  fsorsn; 
pt.indx  =  index_froa_grid(isx+ix, isy+jy , isz-kz-1) ; 
pointCpt_indx3.nl  =  pointCpt_indx3.nl  -  fsorsn; 
pt.indx  *  index_froB_grid(isx+ix,isy-jy,isz+kz) ; 
pointCpt_indx3.nl  «  pointCpt_indx3.nl  +  fsorsn; 
pt.indx  «  index.f roB_grid( isx+ix , isy- jy , isz-kz-1 ) ; 
pointCpt_indx3.nl  *  pointCpt.indx3.nl  -  fsorsn; 
pt.indx  *  index_froB_grid(isx-ix-l, isy+jy, isz+kz) ; 
pointCpt_indx3.nl  •  pointCpt_indx3.nl  +  fsorsn; 
pt.indx  *  index_froB_grid(isx-ix-l, isy+jy, isz-kz-1) ; 
pointCpt_indx3.nl  ■  pointCpt_indx3.nl  -  fsorsn; 
pt.indx  *  index.f roa_grid(isx-ix-l,isy-jy , isz+kz) ; 
pointCpt_indx3.nl  ■  pointCpt_indx3.nl  +  fsorsn; 
pt.indx  *  index.f ron_grid(isx-ix-l,isy-jy, isz-kz-1 ) ; 
pointCpt_indx3.nl  *  pointIpt_indx3.nl  -  fsorsn; 


> 

/*  find  v&lass  at  ths  xy  corners  in  the  z-plsne  of  the  source  */ 

else 

if  ((ix  !*  0)  kk  (jy  !-  0)) 

i 

pt.indx  «  iadex_froa_grid(isx+ix,isy+jy,i*z); 
point Cpt_indx3 .ul  *  point [pt.indx] .ul  +  fsorsn; 
pt.indx  *  index.f ron_grid( isx+ix, isy-jy, isz); 
point [pt.indx] .ul  *  point [pt.indx] .ul  +  fsorsu; 
pt.indx  ■  index_fron_gridCisx-ix, isy+jy, isz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  -  fsorsu; 
pt.indx  «  index_froB_grid(isx-ix,isy-jy ,isz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  -  fsorsu; 

pt.indx  *  index.f ron_grid(isx+ix, isy+jy, isz) ; 
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point [pt.indx] .vl  =  point [pt.indx] .vl  +  fsorsv; 
pt_indx  =  index_froB_grid(i6x+ix,isy-jy-l,isz) ; 
point [pt.indx] .vl  =  point [pt.indx] .vl  -  fsorsv; 
pt_indx  =  index_froB_grid(isx-ix-i,isy+jy ,isz) ; 
point [pt.indx] .vl  =  point [pt.indx] .vl  +  fsorsv; 
pt.indx  =  index_froB_grid(isx-ix-l,isy-jy-l,isz) ; 
point Cpt_indx] .vi  =  point [pt.indx] . vl  -  fsorsv; 

/*  no  u-displaceaent  updates  needed  here  */ 

> 

/*  next  the  xz  comers  uithin  the  y-plane  of  the  source  */ 

else 

if  ((ix  !=  0)  ft*  (kz  !=  0)) 

{ 

pt_indx  =  index_from_grid(isx+ix,isy,isz+kz); 
point [pt_indx] .ul  =  point [pt_indx] .ul  +  fsorsu; 
pt.indx  =  index.f roa_grid(isx+ix, isy ,isz-kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  +  fsorsu; 
pt_indx  =  index_froB_grid(isx-ix,isy,isz+kz); 
point [pt.indx] .ul  =  point [pt.indx] .ul  -  fsorsu; 
pt_indx  =  index.f ron_grid(isx-ix, isy ,isz~kz) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  -  fsorsu; 

pt.indx  *  index.f roa.gr id ( isx+ix, isy, isz+kz ) ; 
point [pt.indx] .ul  *  point [pt.indx] .ul  +  fsorsu; 
pt.indx  =  index_from.grid(isx+ix,isy,isz-kz-l) ; 
point [pt.indx] .ul  =  point [pt _indx] .ul  -  fsorsu; 
pt.indx  =  index.f  ron_gr id ( isx-ix- 1 ,  isy  , isz+kz)  ; 
point [pt.indx] .ul  =  point [pt.indx] .ul  +  fsorsu; 
pt.indx  =  index_fron_grid(isx-ix-l,isy,iaz-kz-l); 
point [pt.indx] .ul  -  point [pt.indx] .ul  -  fsorsu; 

/*  no  v-displaceaent  updates  needed  here  */ 

> 

/*  note:  if  jy  !=  0  and  kz  !=  0  (i.e.  ix  =  0)  then  no  updates 
to  displaceaent  are  needed  since  the  v  and  u  updates 
have  already  been  handled  in  the  previous  sections 
(due  to  the  staggered  grid  geoaetry  and  there  is  no 
u  update  because  of  syaaetry.  */ 

else 

/*  remaining  updates  are  for  the  case  of  tuo  zero  values  and  one 
non-zero  value  for  ix,jy,kz  (but  only  ix  !'  0  is  needed 
because  of  symmetry)  */ 
if  ((ix  ! =  0)  ft*  (jy  *«  0)  ft*  (kz  «  0)) 

pt.indx  *  index.f roa.gr id ( i  ‘ix.isy.isz) ; 
point [pt.indx] .ul  *  pointCp  ,ndx].ui  +  fsorsu; 
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pt.indx  =  index_from_grid(iax-ix,iay ,iaz) ; 
point [pt _ indx] .ul  *  point [pt_indx] .ui  -  f aorau; 


/*  check  on  the  magnitude  of  the  9-3 ,0,0)  forau  tern  - 
if  aba(thia  tara)  <  xoff  then  aource  ia  negligible 
and  ia  turned  off.  lote:  only  check  thia  after  a  full 
cycle  aince  don't  want  to  ahut  off  at  zero  croaaing  */ 
if(ltiae  >  25) 
if (ix  ==  -3) 

if (f aba (f aorau)  <  xoff) 
laora  *  FALSE; 


> 


/*  sssssssssssstsssscMSEcnBttsBsiziiciKiicnBsssismssrESEz 

/* 

/*  Function  atreaa 

/* 

/*  SMIlSSSSZISSSSSSSICtBtlimBtltlUBIIIIISIIIIIIttUlHSSn 

void  aodel: : atreaa (void) 

int  ix,iy,iz; 

ix  «  thia_coordinate(0) ; 
iy  «  thia_coordinate(l) ; 
iz  =  thia_coordinate(2) ; 

if ((ix  <  1  II  ix  >«  (XDIH-1) )  II  (iy  <  1  II  iy  >»  (TOIH-1))  II 
(iz  <  1  II  iz  >■  (ZDIM-1))) 

$ 

elae 

< 

calc.txxdxO ; 
calc_txydy() ; 
calc.txzdz ( ) ; 

calc.txydxO ; 
calc_tyydy(); 
calc.tyzdz ( ) ; 

calc_txzdx() ; 
calc_tyzdy() ; 
calc.tzzdzO; 

> 


*/ 

*/ 

*/ 

*/ 

*/ 
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/* 

/*  Function  update 

/* 

/*  s:ssss:sssssscsssssscssssssss:sssssssss=sss: 

void  nodel: :update(void) 

u2  =  ul; 
v2  *  vi; 
v2  =  vl; 

ul  =  uO; 
vl  =  vO; 
ul  >  lO; 

> 

/* 

/* 

/* 


End  of  vave_3d  functions 


/*  Bava.3d.h8  -  include  iile  for  Bava_3d.cs  */ 

/*  */ 

/*  Programmer:  Julia  Allan  */ 

/*  Data:  28  February  1989  */ 

/*  */ 

/*  */ 

/*  Copyright  (c)  1989  */ 

/*  Information  Processing  and  Communications  Laboratory  */ 

/a  Voods  Hole  Oceanographic  Institution  a/ 

/*  All  rights  reserved.  */ 

/*  */ 

/a  This  material  cannot  be  distributed  or  sold  Bithout  a/ 

/a  prior  permission  of  the  author(s).  a/ 

/*  */ 


/*  :nn»:inun»tss<»itsi»i»»HE>ict»nt»ncEinnnst  * / 


/•  define  macros  for  timing  a/ 

•define  TINER.OI  CM.start.timerCl) 

•define  TINER.RESET  CN.reset.timerO 

•define  TIHER.OFF  (void)  CM.stop.timer(l) 


/a  common  definitions  for  3D  finite  difference  modeling 


•define 

ZDIN 

•define 

YDIN 

•define 

ZDIN 

•define 

TOTAL 

•define 

NTEL 

•define 

DEL 

•define 

PI 

•define 

TRUE 

•define 

FALSE 

•define 

SET.DEBUG 

64  /a  model  dimensions 

64 

64 

ZDINaYDIlfaZDIM 

4 

6 

3. 14159266f 
1 
0 

FALSE 


*/ 

*/ 


/a  define  common  variables  */ 

mono  char  out.fb  *  FALSE,  out.xy  ■  FALSE,  out.yz  *  FALSE,  out.xz  *  FALSE; 

mono  int  snap_file,snap_fb,nt,ntl; 

mono  float  zlamb,ylamb,zlamb,zylamb,zzlamb,yzlamb; 

mono  float  dz.dy ,dz,dt,fpeak; 

mono  int  isz.isy.isz; 

mono  char  lsors,  DEBUG  *  FALSE; 


/a  define  the  domain  for 
domain  model  { 

float  u0,ul,u2; 
float  v0,vl,v2; 
float  b0,b1,b2; 


the  parallel  variables  -  27 


variables  =  108  bytes  a/ 
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float  clanb; 
float  can; 
float  clp2au; 
float  rho; 

float  txxdx.txydy ,txzdz; 
float  txydx.tyydy .tyzdz; 
float  txzdx.tyzdy.tzzdz; 
float  alpha; 

float  div.curlx, curly, cur lz; 

float  alpha.f ct(int  indx) ; 

void  boundary (void) ; 

void  calc_adiap(float  *adisp) ; 

void  calc.bdiap (float  *bdisp) ; 

void  calc_cdiap(float  *cdi«p) ; 

void  calc_displ(void) ; 

void  calc.txxdx(void) ; 

void  calc_txydy(void) ; 

void  calc.txzdz(void) ; 

void  calc.txydx(void) ; 

void  calc_tyydy(void) ; 

void  calc.tyzdz(void) ; 

void  calc.txzdx(void) ; 

void  calc.tyzdy(void) ; 

void  calc.tzzdz(void) ; 

void  divcurl(void) ; 

void  initialize (void) ; 

void  output_fb(void) ; 

void  atreaa(void) ; 

void  updata(void) ; 

>  point [TOTAL] ; 
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D.2  convert  program 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c===== 

c 

c 

c 

c 

c 

c 


(SSSU 


c 


c 

1001 

c 


c 

c 

c 

c 

c 

c 

c 


program  convert 


Programmer:  Julie  Allen 
Date:  16  June  1989 


Copyright  (c)  1989 

Information  Processing  and  Communications  Laboratory 
Hoods  Hole  Oceanographic  Institution 
All  rights  reserved. 

This  material  cannot  be  distributed  or  sold  without 
prior  permission  of  the  author(s). 


»:::<»tts:::nr»tnniu»«i»>iaim>isnnn>tRK:Kt 

Program  is  currently  hardwired  for  models  with  dimensions 
64x64x64 

Convert  formatted  files  output  on  CM  and  transferred  via  ftp  to 
the  VAX  to  binary  files  ready  for  input  to  program  SIAP 


real *4  value (64 ,64) 
integer*4  inf ile , ioutf ile 

data  inf ile/11/, ioutf ile/ 12/ 

format (a) 

print  '(lx,a,$) 'Enter  scale  factor  (1  for  no  scaling):  ' 
read  *, scale 

Open  input  file 

open ( inf ile , status* ’ old ’ , form* ’ f  ormatted ' ) 

Open  output  file 

open (ioutf ile .status* ’new ' .form* 'unformatted' ) 

read (inf ile, 100 1 ,err*800,end*900)  dummy 
do  1*1,64 

do  j  *  1,64,8 

read(infile,*,err*800,end*900)  ( valued, k)  ,k*j  ,  j+7) 
enddo 
enddo 
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c 

c 

c 


c 


c 

c 

800 

c 

900 


c= 


Scale  the  data 

if (scale. ne. 1.0)  then 
valmin  *  1000.0 
v aim ax  =  -1000.0 
do  i=l,64 

do  j=l,64 

valued, j)  =  valued, j)  *  scale 
if  (valued,  j)  .It  .valmin)  valmin  =  valued,  j) 
if  (valued,  j)  .gt.valmax)  valmax  =  valued,  j) 
enddo 
enddo 

print  *a 'minimum:  ’, valmin,*  maximum:  ’, valmax 
endif 

write (ioutf ile) 
write (ioutfile) 

write(ioutfile)  ((valued,  j)  ,i=l  ,64)  ,  j=l  ,64) 
go  to  900 

call  err sns ( i err , , , inn it ) 

print  * , *  errsns  #  ’ , ierr , *  on  unit  ’ , iunit 

close (inf ile) 

close(ioutfile) 

stop 

end 

:::sss:sts3sscs«9ess«sssS8ssssssa*sc8SSi»«B8scssessssssx 
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D.3  n-dimensional  grid  library 


The  following  source  code,  used  in  the  wave_3d  program,  was  written  by  Robert  Whaley,  Site 
Representative  from  TMC  at  the  NRL  Connection  Machine  Facility. 
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void  void:  :uii_|rid(MBo  lot  z,  Mono  int  y,  aono  int  c,  mono  int  a,  int  void: :  next  (mono  int  dia,  int  void::  *  aono  value.to.get)  { 

aono  int  b,  mono  int  c,  mono  int  d)  {  int  return. value; 

aono  unaigned  diaeneion_array[7] ;  aono  CK.aeaaddr.t  kludge  -  value. to.gat ; 

diaeneion.arrayCO]  •  z;  CM.get.froa.nese.lLOretum.value,  kludge,  dia,  CM.upsard, 

dlaenaion. array  [1]  *  y;  bit.eizeof (retum.value))  ; 
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domain  void  •  *oid: : pointor_from.gr id(int  z)  { 

\*  nd_grid . ha  -  inclod*  ill*  for  n-dia*n.ional  grid  library  */  r«tnm(CUBEiDDR_TO_POI»TER(domain  »oid,  ind«x_from_grid(x))) ; 
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/*  fb  -  C*  frame 

buffer  functions 

*/ 

/* 

*/ 

/*  Author: 

Robert  Whaley 

*/ 

/*  E-nail: 

ehaleyCnrl-caf . arpa 

*/ 

/*  US  nail: 

Robert  Whaley 

*/ 

/* 

Vaval  Research  Lab 

*/ 

/* 

Code  5150 

*/ 

/* 

Washington,  O.C.  20376-S000 

*/ 

/*  Phone: 

(202)  404-7019 

*/ 

/*  Affiliate: 

TMC  (Site  Rep  at  IRL) 

*/ 

•include  <stdio.hs> 

•include  <cn/paris .hs> 

•include  <c«/cmfb.hs> 

•include  "nd.grid.hs" 

•include  "fb.hs" 

•define  bit.sizeof (x)  (8  *  sizeof(x)) 

static  struct  CMFB.display.id  ay .display; 
static  CK_vp_set_id_t  ny.vp.set ; 
static  CM_vp_set_id_t  old.vp.set; 
static  CM_vp_sst_id_t  sane_geo.vp.set; 
static  CM_geo*etry_id_t  ny_geo«etry ; 
static  CM.aeaaddr.t  ny .color; 

int  CMFB_width( struct  CMFB_display_id  *); 
int  CKFB_height( struct  CMFB.display.id  *); 

CHFB.buff er_ id.t  CMFB.spare.buff er( struct  CMFB.display.id  *) ; 
void  CHFB_svitch_buffer( struct  CMFB.display.id  *,  CM.buff er.id.t) ; 
unsigned  .CMI.get.f ield.address.fron.f ield_id_safely(CM_f ield_id_t) ; 
unsigned  .CMI.f ield_location(unsigned) ; 

void  .CMI.physical.nove.aleays (unsigned,  unsigned,  unsigned); 
char  *getenv(char  *  nane) ; 

void  init_frane_buffer(int  x_dim,  int  y.din)  { 
unsigned  dins [2] ; 
int  physical.!,  physical.y; 
int  zoon; 

■ane_geo_vp_set  =  CM.allocate.vp.set (CM. vp.set.geosetry (CM_ current. vp.set) ) ; 
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dims [0]  =  x.dim; 
dims  Cl]  =  y.dim; 

my.geometry  =  CM_create_geometry(dims,2) ; 
my_vp_set  =  CM.allocate.vp.set (my.geometry) ; 

if  (CMFB.attach.display (getenv("CM_FRAMEBUFFER") ,  toy .display) )  { 
if  (getenv("CM_FRAMEBUFFER") )  { 
printf  ("imt.frame.buffer:  Could  not  attach  to  frame  buffer: 

*/.s.\n", 

getenv("CM_FRAMEBUFFER")) ; 

printf ("Framebuffer  is  probably  not  connected  to  your  sequencers  An") ; 

> 

else  { 

printf ("init.frame.buffer :  Could  not  attach  frame  buffer  An"); 
printf ( "Probably  using  sequencers  that  have  no  frame  buffer An"); 

> 

> 

CMFB_initialize_display(tay .display ,  8,  1); 

CMFB-.  ini  t  ial  ize.  color  .table  ( toy  _di  splay  ) ; 
physical.!  =  CMFB_vidth(tay .display) ; 
physical.y  =  CMFB.height (toy .display ) ; 

zoom  =  ((physical.*  /  x.dim)  <?  (physical.y  /  y.dim))  -  1; 
CMFB_set_zoom(tay_display ,  zoom,  zoom,  0); 
old.vp.set  *  CM.current.vp.set ; 

CM.aet.vp.set (my.vp.set) ; 

*y_color  »  ( CM.memaddr _t )  CM.allocate.heap.f ield(8) ; 
CM_u_nove_zero_always_lL(my_color ,  8) ; 

CM. a et.vp. set (old.vp.set) ; 


void  release.frame.buff er()  { 

CMFB_detach_display(tmy .display) ; 
CM.deallocate.heap.f ield (my. color) ; 
CM.deallocate.vp.set (my.vp.set) ; 
CM.deallocate.geometry (my .geometry) ; 


void  void: :plot_x_y (unsigned  short  x,  unsigned  short  y, 

unsigned  char  color)  { 
mono  CM.vp_set.id_t  old.vp.set; 
mono  CMFB.buf f er.id.t  the.buffer; 
unsigned  int  a.send.address ; 
old.vp.set  *  CM.current.vp.set; 

CM. set.vp. set (my.vp.set) ; 
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CM_u_move_zero_alHays_lL(my_color ,  8) ; 

CM_set_vp_set(old_vp_set) ; 

CMFB_shuffle_from_x_y(ka_send_address,  kx,  ky,  my .geometry) ; 
CM_send_lL(my_color ,  fca_send_address ,  kcolor,  bit.sizeof (color) , 

( CM_f i eld_ id_t ) CM_do_not  _not if y .token ()) ; 

CM.set.vp.setC  ny.vp.set) ; 

the.buffer  =  CMFB.spare.buf f er(kmy.display) ; 

CMFB_urite_preshuffled_alHays(kmy_display ,  the.buffer,  my.color, 

0,  0); 

CMFB.snitch.buff er (kmy.display ,  the.buffer) ; 

CM_set_vp_set(old_vp_set) ; 

> 

void  void: :plot_x_y_over (unsigned  short  x,  unsigned  short  y, 

unsigned  char  color)  { 
mono  CM.vp_set.id_t  old.vp.set; 
mono  CMFB.buffer.id_t  the.buffer; 
unsigned  int  a.send.address ; 
old.vp.set  =  CM.current.vp.set ; 

CM_set_vp_set(my_vp_set) ; 

CM_set_vp_8et(old_vp_set) ; 

CMFB_shuffle_from_x_y(ka_send_address,  kx,  ky,  my .geometry) ; 
CM_send_lL(my_color,  ka.send.address ,  kcolor,  bit.sizeof (color) , 

(CM.f ield. id.t ) CM.do.not.not if y .token ()) ; 

CM.set.vp.set (my.vp.set)  ; 

the.buffer  =  CMFB.spare.buf fer (kmy.display) ; 

CMFB.Hrite.preshuffled.alHays (kmy.display ,  the.buffer,  my.color, 

0.  0); 

CMFB.ssitch.buf fer (kmy.display,  the.buffer) ; 

CM_set.vp.set (old.vp.set) ; 

> 

void  void: :plot_fro*_grid(unsigned  char  color)  { 
mono  CMFB.buf fer. id.t  the.buffer; 
mono  int  vp_ index; 
void  void::  *  mono  color.prime; 

mono  unsigned  color .phys.loc,  color _prime.phys.loc; 
color.phys.loc  *  (unsigned)  kcolor; 
old.vp.set  *  CM.current.vp.set; 

CM_set_vp_set(same_geo_vp_set) ; 

vp.index  ■  CM.geometry .total. vp.ratio (CM. vp.set.geometry (old.vp.set)) ; 
color.prime  =  CM.allocate.stack.f ield(8) ; 
color _prime.phys.loc  = 

.CMI.f ield_location( 
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.CMI.get.l ield.addres s.lrom.l ield.id.sal ely ( color .prime) ) ; 
while  (vp.index —  >  0)  •{ 

.CMI.physical.move.always (color _prime.phys.loc,  color _phys_loc, 

8); 

color.prime.phys.loc  +=  8; 

color_phys_loc  +=  (unsigned)  CM.user.memory.address.limit  +  4; 

> 

the.buffer  =  CMFB_spare_builer (lnny_display) ; 

CMFB.write.always (toy .display ,  the.buller,  color.prime,  0,  0); 
CMFB_switch_buffer(toy .display ,  the.buller) ; 
CM_deallocate_stack_through( color .prime) ; 
CM_set_vp_set(old_vp_set) ; 


/*  Ib.hs  -  include  tile  lor  lb  library  */ 
void  set_color(int  color.id,  int  red,  int  green,  int  blue)  { 
CMFB.write.color (tay_di8play ,  CMFB.red,  color.id,  red); 
CMFB.write.color (toy .display,  CHFB_green,  color.id,  green); 
CHFB_write_color( toy .display,  CMFB.blue ,  color  id,  blue); 

> 

void  init_lrame_buller(int  x.dim,  int  y.dim) ; 
void  release_lrame_buller(void) ; 

void  void: :plot_x_y (unsigned  short  z,  unsigned  short  y, 

unsigned  char  color); 

void  void: :plot_x_y .over (unsigned  short  x,  unsigned  short  y, 

unsigned  char  color); 

void  void: :plot_lro«_grid(unsigned  char  color); 

void  set_eolor(int  color.id,  int  red,  int  green,  int  blue); 
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